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We present quantum algorithms for the simulation of quantum systems in one spatial dimension, 
which result in quantum speedups that range from superpolynomial to polynomial. We first describe 
a method to simulate the evolution of the quantum harmonic oscillator (QHO) based on a refined 
analysis of the Trotter-Suzuki formula that exploits the Lie algebra structure. For total evolution 
time t and precision e > 0 , the complexity of our method is 0 (exp( 7 -y/log(A/e))), where 7 > 0 is 
a constant and N is the quantum number associated with an “energy cutoff” of the initial state. 
Remarkably, this complexity is subpolynomial in N/e. We also provide a method to prepare discrete 
versions of the eigenstates of the QHO of complexity polynomial in log{N)/e, where N is the 
dimension or number of points in the discretization. This method may be of independent interest as 
it provides a way to prepare, e.g., quantum states with Gaussian-like amplitudes. Next, we consider 
a system with a quartic potential. Our numerical simulations suggest a method for simulating the 
evolution of sublinear complexity for constant t and e. We also analyze complex 

one-dimensional systems and prove a complexity bound 0{N), under fairly general assumptions. 
Our quantum algorithms may find applications in other problems. As an example, we discuss 
the fractional Fourier transform, a generalization of the Fourier transform that is useful for signal 
analysis and can be formulated in terms of the evolution of the QHO. 

PACS numbers: 


I. INTRODUCTION 

One of the best known problems that a quantum com¬ 
puter is expected to tackle more efficiently than a classi¬ 
cal one is quantum simulation (QS) [TH3j. QS provides in¬ 
sight to the behavior of a quantum system, hence it serves 
as a powerful method applicable to a variety of scien¬ 
tific areas including physics (c.f., [HE]), quantum chem¬ 
istry (c.f. |6H9]), and computer science (c.f., [IQ]). One 
way to explore this problem is by solving Schrddinger’s 
equation, so that the state of the system evolves under 
the evolution operator, as determined by the Hamiltonian 
of the system. A quantum simulator mm carries out 
this simulation task on a quantum system or computer 
by an approximate implementation of the evolution, ei¬ 
ther as a sequence of elementary gates (i.e., a digitial 
quantum simulator) or by direct implementation of the 
Hamiltonians (i.e., an analog quantum simulator); this 
paper is concerned with digital quantum simulators. 

Simulating a quantum system requires a certain 
amount of “time and memory”, called resources, and the 
important question is how the resources scale with quan¬ 
tities such as the volume or size of the system (or the di¬ 
mension of the Hilbert space), and precision. Generally, 
when the simulation is carried on a classical computer, it 
would require an undesirably large amount of resources, 
which R.P. Feynman addressed as an “exponential explo¬ 
sion” , and it rapidly becomes an intractable problem even 
for relatively small systems. On the contrary, a quantum 
simulator is expected to carry this simulation efficiently, 
at least in some important cases [T]. 

A QS may be implemented for two goals: i- to com¬ 
pute time-dependent physical properties, such as scat¬ 


tering amplitudes, as determined by the evolution of the 
quantum system, and ii- to compute spectral properties, 
such as eigenvalues or expectation values of observables 
on different eigenstates of the system. In the first goal, 
we are mainly concerned with the simulation of the evo¬ 
lution operator U{t) := exp{—iHt}. Results on the so- 
called Hamiltonian simulation problem provide efficient 
quantum algorithms to simulate U{t) on quantum com¬ 
puters da [MU]. Generally, the complexity of such al¬ 
gorithms (i.e., the number of elementary two-qubit gates 
and queries) depends only polynomially in r, d, and 1/e, 
where t = ||i7t||, d is the sparseness of H, and e > 0 
is a precision parameter. The Trotter-Suzuki approxi¬ 
mation [141119] plays an important role in these results. 
The complexity dependence on 1/e can be exponentially 
improved using recent methods that implement a series 
approximation ofU{t) [SOlEI]. In the second goal, we are 
mainly concerned with the preparation of eigenstates iia- 
125] . Two commonly used techniques for these quantum 
algorithms are the simulation of quantum adiabatic evo¬ 
lutions |2^ and the implementation of the so-called phase 
estimation algorithm m- But since approximating the 
ground state energy of many-body systems is a computa¬ 
tionally hard problem |28| . quantum algorithms to pre¬ 
pare ground states are generally inefficient. 

The quest for more efficient methods for QS is ongo¬ 
ing, particularly sparked for the need to understand the 
quantum advantages of relatively small, potentially near¬ 
future, quantum computers (c.f., [22]). Unfortunately, 
most quantum algorithms for QS up to date concern the 
case of discrete, finite dimensional quantum many-body 
systems where ||iL|j < 00 . These algorithms cannot be 
directly applied to speedup the simulation of continuous- 
variables (GVs) quantum systems: the effective energy 
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scale of the problem can increase polynomially with the 
dimension of the relevant Hilbert space. While a few ex¬ 
ceptions exist (e.g., [551130]), a detailed analysis of the 
power and limits of quantum computers for QS of CV 
systems is lacking. Here, we incentivize such an analysis 
by studying a quantum-computer simulation of simple 
one-dimensional quantum systems. 

We first consider the quantum harmonic oscillator 
(QHO). Remarkably, our first result is a quantum al¬ 
gorithm that simulates U(t) with subexponential com¬ 
plexity. That is, the complexity is 0(exp(7-\/log(A^/e))), 
where N is the dimension (as determined by the dis¬ 
cretization and effective energy of the problem) and e 
is the precision. This complexity does not depend on 
the evolution time because the evolution operator of the 
QHO is periodic. The first result is obtained via a re¬ 
fined analysis of the Trotter-Suzuki approximation that 
exploits the Lie algebra structure of this problem, allow¬ 
ing us to bypass the no fast-forwarding theorem of m, 
which would indicate a complexity that is at least linear 
in N. Our results are in contrast with, for example, the 
results in [551 130] i where the complexity is polynomial 
in N. Our result suggests a quantum superpolynomial 
speedup over the corresponding classical algorithms for 
this problem, showing a significant quantum advantage 
even for the simulation of this simple model. 

We then present quantum algorithms of complexity 
polynomial in log(A'^)/e that prepare approximations of 
the eigenstates of a discrete version of the QHO. These 
states are prepared by simulating the evolution opera¬ 
tor of a Hamiltonian that is a discrete version of the 
Jaynes-Cummings model, using the Trotter-Suzuki ap¬ 
proximation. In particular, for the preparation of the 
eigenstate of lowest eigenvalue, the complexity is polyno¬ 
mial in log(A^/e). The computation of spectral properties 
of the QHO can then be done by computing expectation 
values on such states using a variety of known techniques 
(c.f., [31] )■ These algorithms may be of independent in¬ 
terest as they allow for a very efficient way to prepare 
states whose amplitudes are, for example, Gaussian-like. 

In order to understand the complexity of simulat¬ 
ing more complex quantum systems, we perform sev¬ 
eral numerical simulations of a quantum system with 
a quartic potential. In contrast with the QHO, quan¬ 
tum algorithms to simulate U{t) based on high-order 
Trotter-Suzuki approximations seem to be of complex¬ 
ity sublinear in N for this case. In particular, our 
numerical simulations suggest that such a complexity 
is for arbitrarily small, but con¬ 

stant, 77 > 0. The O notation hides factors that are 
polynomial in log(A^|t|/e). This complexity represents a 
polynomial quantum speedup, with respect to N, over 
the classical algorithm. The sublinear complexity in N 
may be a result of the algebraic structure satisfied by 
the operators in the Hamiltonian (see the recent results 
in [35] for more details). 

Finally, we use our results in m to present a generic 
quantum method to simulate the dynamics of one dimen¬ 


sional quantum systems of complexity 0(|t|7V), under 
fairly general assumptions. In contrast with our previ¬ 
ous results, the method in [5T] implements a series ap¬ 
proximation of the evolution operator and the complexity 
dependence on 1/e is polylogarithmic. Our main contri¬ 
bution in this case is to represent the evolution operator 
in the interaction picture, so that potentials of high de¬ 
gree do not increase the complexity significantly. We note 
that classical algorithms for this problem are expected to 
have complexity that is super-linear in N. 

The remainder of the paper is organized as follows. In 
Sec. [IT] we revisit the QHO where we discuss particular 
properties that are useful for our quantum algorithms, 
and also classical algorithms for simulating this model. In 
Sec. |HI| we define a discretization of the QHO as the start¬ 
ing point for a QS. We provide a quantum algorithm to 
compute scattering amplitudes associated with the QHO 
in Sec. jlVj and quantum algorithms to prepare the eigen¬ 
states of the discrete QHO in Sec. [V] We analyze the 
quantum system with a quartic potential in Sec. jVIj and 
present the upper bound on the complexity of simulating 
general one-dimensional quantum systems in Sec. [VH] In 
Sec. [VHll we discuss related work with particular empha¬ 
sis on the so-called fractional Fourier transform, which is 
a generalization of the Fourier transform used in signal 
analysis [33]. We state the conclusions in Sec.]TX 


II. THE QHO REVISITED 

The QHO is described by its Hamiltonian {h = 1) [34] . 

H=^[x^ + {-zd,r]. ( 1 ) 

Because x G (— 00 , 00 ), we sometimes refer to H as the 
continuous-variable or CV QHO. The eigenfunctions of 
H are the so-called Hermite functions, 

1 

■Ipnix) = . e 2 llr,{x) , (2) 

where each Hn{x) is the n-th (physicists’) Hermite poly¬ 
nomial and n = 0, 1, 2,... are quantum numbers. Then, 

Hipn{x) = {n+ l/2)'0„(x) , (3) 

where n-|-l/2 are the corresponding eigenvalues. In stan¬ 
dard bra-ket notation, we represent the eigenfunctions 
ipn{x) by the eigenstates |'(/„) and the QHO, in operator 
form, is 

H=l{x^+p^). (4) 

Then, iLlV'n) = {n + l/2)\ipn), and the position and mo¬ 
mentum operators satisfy 


{ll^m\x\'tpn) = 

j dx 'tljrn{x){xipri{x)) , 

( 5 ) 

{lpm\p\tpn) = 

-i J dx 1pm{x){da;^n{x)) , 

( 6 ) 
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respectively. The evolution operator of the QHO for time 
t is U{t) = exp{—iHt}. Since the evolution operator is 
periodic, i.e. C/(0) = C/(47r), we can assume < = 0(1)- We 
also use \x) to denote the eigenstates of x of eigenvalue 
X G (— 00 , 00 ). 

We can alternatively define the raising and lowering 
operators, 


in terms of IV'n), to obtain the amplitudes c„ and c'„. 
Because the sum in Eq. ([^ involves 0{N') terms, the 
worst-case complexity of the classical method is of order 
polynomial in iV' if < 00 . We also note that 

= J J dx'dx 'tpm{x')A{x',x)lpn{x) , (9) 


i' = {x — ip )/, a = {x + ip)/V 2 , 


(7) 


and write H = a^a-|-l/2. The eigenstates satisfy 
a^a\ipn) = and a^\\j}r/) = Vn + l|^/>„+i). Then, 

other eigenstates of the QHO can be obtained by re¬ 
peated action of on the vacuum (ground) state |'0o)- 
Later, we will use this property to devise a quantum al¬ 
gorithm that prepares eigenstates, up to some approxi¬ 
mation error (Sec. |V). 

A well-known resu t states that x and p are related via 
the Fourier transform (FT), which transforms x —>■ —p 
and p ^ X (i.e., a 7r/2 rotation in phase space). Simi¬ 
larly, the FT transforms —>■ —ia^ and a —>■ ia. One 
implication is that the eigenstates \ijjn) of the QHO are 
also eigenstates of the FT, and direct computation shows 
that the eigenvalues of the FT are (—*)". We will use 
this property to prove some results regarding the quan¬ 
tum algorithm that simulates the evolution of the QHO 
(Appx. A). 

It is also important to remark that the operators in 
H generate a Lie algebra sp{2) of dimension 3. This 
property will be useful to bound the errors when approx¬ 
imating the evolution operator in Appx. B. The result is 
that the effective norm of nested commutators is signif¬ 
icantly smaller than the product of the effective norms, 
allowing us to perform a refined analysis of the errors in 
Hamiltonian simulation methods. 

Of particular interest in a QS is the computation of 
quantities such as scattering amplitudes, {(p'\ip{t)). Here, 
|(p(<)) := U{t)\(p) is the evolved state and \‘p), \‘p') are 
some other specific states of the system, such as eigen¬ 
states of the position or momentum operator, or more 
general states. Also important is the computation of 
expectation values or correlation functions in different 
eigenstates, namely (■i/'ml A|i/'n), where A is some observ¬ 
able; e.g., A = p’-^x^^, where li,l 2 > 0. Our quantum 
algorithms are designed to compute such quantities, but 
they could also be used for the computation of more gen¬ 
eral quantities, such as expectation values of other uni¬ 
tary operators, after minor and straightforward modifi¬ 
cations. 

With no loss of generality, \p) = Y/n=o^'n-\'>Pn) and 


W)=T.n=,<^'n\/^n) (with EJCnP=E 
that 


„ iC'nP = 1), SO 


N' 


{p'\p{t)) = 


n=0 


( 8 ) 


where A{x’,x) = {x'\ A |a;). Assuming n,m < N', classi¬ 
cal methods to approximate the integral in Eq. can 
also have a worst-case complexity of order polynomial in 
N'. This is because the functions 'ijjn{x) present oscilla¬ 
tions that become more significant as n increases, so a 
good approximation to the integral by a finite sum can 
only be obtained if the number of terms in the sum is 
polynomial in N' [2H]. Nevertheless, it is well-known that 
many quantities associated with the QHO can be analyt¬ 
ically obtained, so our quantum algorithms will be more 
powerful than classical ones only in certain scenarios. 

A natural quantum method to compute Eq. ([^ in¬ 
volves direct QS of the QHO. In this case, if the states 
\p) and \lp') can be efficiently prepared, the quantum 
method is efficient or not whether U(t) can be simulated 
efficiently or not, respectively. In this paper we are first 
interested in devising efficient quantum-computer sim¬ 
ulations of the evolution induced by the QHO and then 
we will consider quantum algorithms for preparing eigen¬ 
states. 


III. A DISCRETE QUANTUM HARMONIC 
OSCILLATOR 


In analogy with the CV QHO, we define a discrete 
QHO by the Hamiltonian 

77 " + (/)")• ( 10 ) 


The Hilbert space dimension is N, where A > 2 is even 
for simplicity, x^ is the discrete “position” operator given 
by the N x N diagonal matrix 


x'^ = 


(-N 


0 


27r 1 

A 2 


0 i-N+ 2) 


V 0 


0 


(TV-2)7 


( 11 ) 


and p‘^ is the discrete “momentum” operator given by 


p^ = {F^)-\x‘^.F^ . ( 12 ) 


The N X N unitary matrix F// is the so-called centered 
discrete Fourier transform, which is the standard discrete 
Fourier transform F*^ up to a simple (cyclic) permutation. 
Its matrix entries are 


Then, a standard classical method to compute scattering 
amplitudes involves the spectral decomposition of \ip) and 


fe = 27rjA:/A) , (13) 
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where j,k G {—JV/2,...,JV/2 — 1} label the rows and 

columns, respectively. Then, 

and 


X = 


/O 1 0 
'O 0 1 


0 0 0 

Vl 0 0 




1 

o/ 


(14) 


is the operation that performs a cyclic permutation, shift¬ 
ing the indices by one. The relation between and 
will be useful to prov ide an efficient quantum circuit that 
implements F^'^ (Sec. IVA). 


We will generally assume that we are in the limit of 
large N, although some results do still apply when N is 
fairly small. Then, in the following, the order notation 
to bound approximation errors assume the asymptotic 
limit; we refer to the corresponding appendices for more 
details. 


A. Spectral properties 


In contrast to the CV QHO, the Hamiltonian may 
not be exactly solvable. However, some spectral proper¬ 
ties of can be well approximated from those of the CV 
version. We write and |<()([) for the eigenvalues and 
eigenstates of respectively, and n = 0,1,...,V — 1. 
We also introduce the (unnormalized) quantum states 



which, as we will show, approximate the eigenstates of 
Here, Xj = j^2TT/N and tpn{x) is the n-th Her- 
mite function, so that represents a “discrete Hermite 
state”. In Appx. A, Cor. we show that there exists a 
constant c, 1 > c > 0, such that 

-{n + I/2)|V^d)f = exp(-H(iV)) , (16) 


for all n < cN. Here, for a matrix A, |jA|| is the spectral 
norm and || |^) H is the Euclidean norm of a state |^) . 
The notation exp(—H(V)) states that there is a constant 
/3 > 0 such that the right hand side of Eq. (16) is at most 
g-/3Ar, sufficiently large N. 

Since c < 1, we refer to the subspace spanned by !?/’()), 
for all n < cN, as the “low-energy” subspace. Intuitively, 
in such a low-energy sector, the Hermite functions can 
be well approximated by piecewise constant functions 
'4’n{x) = ipn{xj) if Xj < X < Xj + 2ttJN. The inte¬ 
grals needed to compute properties of the QHO can be 
replaced, within high accuracy, by the sums that appear 
in the discrete case. In contrast, for large values of n, 
the Hermite functions may present oscillations that will 
not be captured under such an approximation, and the 
approximation error gets large. 


The states |'05)) form almost an orthonormal basis of 
the low-energy subspace. In Appx. A, Cor. we prove 

- ^n,m\ = exp(-H(V)) , (17) 

for all n,m < cN. Equations (161 and 0 imply that 
the eigenvalues of in the low-energy sector satisfy 

|Fd-(n +1/2)1 =exp(-H(V)). (18) 

This follows from noticing that {H'^ — {n + 1/2))^ is non¬ 
negative and its smallest eigenvalue is bounded from be¬ 
low by 0 and bounded from above by — (n + 

1/2))^|'0(()/|||'!/(J)P = exp(—H(V)). The corresponding 
eigenstates of satisfy 


|||</(l)-|V>)))||=exp(-H(Ar)) 


(19) 


when n < cN. Since the eigenvalues in the 
low energy subspace are gapped, Eq. (191 can 
be shown by considering the projector = 

lim„_>oo . In particular, for precision 

exp(—0(V)), it suffices to choose u = 0{N). In this 
case, = (I — ds — 

and then “ 1| = 

exp(-H(V)) or ||||0(</^|iA^)|| - 1| = exp(-H(V)). 

The values of the constants hidden by the order no¬ 
tation, as well as c, may be estimated from bounds on 
the tails of the Hermite functions and the correspond¬ 
ing proofs in Appx. A. The constant in the exponen¬ 
tials of Eqs. 0, 0, and ( [I^ may depend on n (i.e., 
it decreases as n increases) but, for our results, it suf¬ 
fices to claim that there is a constant such that the 
approximations of the eigenvalues and eigenvectors are 
exp(—r2(V)). Additionally, Lemmarequires c < 7r/16 
and the corresponding constant in the exponential for 
this lemma is 7r/2 when n = 0. Nevertheless, we can 
also perform a simple numerical analysis to validate our 
results and give better estimates of such constants and 
c. As an example, in Fig. we show the absolute value 
of the overlap between and the actual eigenstates of 
that is, |, for n, m G {0,1,..., N—1}. These 

numerical computations suggest that the states are 
excellent approximations of the eigenstates of for a 
fraction c > 3/4 of the whole spectrum. 

In Fig. [^we show the eigenvalues of the CV QHO and 
H^. Our numerical results suggest that |F() — (n + 1/2)| 
approaches zero as long as n < cN, also for some c > 3/4. 


In Fig. we show the absolute difference between the 
n-th eigenvalue of and n + 1/2, for n = N/2 and n = 
3V/4. The numerical results indicate that this difference 
decays exponentially with N, as determined in Eq. (181, 
and where the constant in the exponential of Eq. (18) 
depends on n. 


IV. SCATTERING AMPLITUDES 

We address the first goal of a QS. In particular, we 
seek a quantum algorithm, built upon a sequence of 
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FIG. 1: Overlap between the eigenstates of a-nd 

the discrete Hermite states, \ipn), for dimension N = 800. 



FIG. 2: Eigenvalues of H’^ for N = 400 (purple) and N = 800 
(yellow), and comparison with n+ 1/2, i.e. the eigenvalues of 
the GV QHO (blue). 


two-qubit gates, to compute scattering amplitudes as 
determined by the evolution of the CV QHO. We let 
:= exp{—iH'^t} be the unitary evolution opera¬ 
tor of the discrete QHO. Intuitively, U^{t) approximates 
U{t), in some sense, as the dimension N grows larger. 
Our goal is to compute {(p'\U(t)\ip), where we assume 
Iv^) = < oo- We also assume that 

W) = TZ=oc'Mn) or \ip') = {2-K/NY/^\xj), where \xj) 
is the eigenstate of the position operator x with eigen¬ 
value . (Note that the eigenstates of x can¬ 

not be written as a finite linear combination of I'f/’n)-) 
For the discrete case, we define oc J2n=o 

and oc J2n=ocMn) or = |j), depending on 

the case. With these definitions, {xj\(p) oc and 

{xj\ (p') oc {j\ so that the states and cor¬ 
respond to discretizations of \p) and \p') in space, re¬ 
spectively. The proportionality constants are needed for 
normalization. The first result of this section is: 

Theorem 1. Let t be the evolution time and e > 0 &e 


(a) 




FIG. 3: Plot of the natural logarithm of the absolute dif¬ 
ference between the n-th eigenvalue of and n -|- 1/2, as a 
function of the dimension (red dots). The linear fit (blue line) 
indicates that such a difference decays exponentially with N. 
(a) n = N/2. The approximate constant in the exponential 
of Eq. (181 is 0.248 in this case, (b) n = 3N/ 4. The approx¬ 
imate constant in the exponential of Eq. (181 is 0.010 in this 


a precision parameter. Assume |t| > 1 and t = 0(1) 
with no loss of generality. Then, there exists N = 
0(log(l/e) -\- N’) such that 

\{p'‘^\U‘^it)\p^)-{p'\U{t)\p)\=0{e). (20) 

Proof. The result follows from the subadditivity property 
of errors. First, we assume oc J2n=o 
a and of be the constants of proportionality, so that 


N' 

n,n'—0 

( 21 ) 

Because T,n=o = T,n=o \^'n? = 1. Eq. @ implies 
|q! — 1| = exp(—H(Af)) and \a' — 1| = exp(— H(7V)). This 
result assumes a choice for the dimension of the discrete 
system of iV = 0{N'), so that 0 < n < N' < cN. The 
value of n within the range represents the low-energy 
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subspace, and we also used that J2n=o “ 0{\/N), 

J2n=o\^n\ = 0{VN). Furthermore, Eqs. ( [T7| and 
imply 

= \t\exp{-n{N)) . (22) 

Note that we can assume t = 0(1) since U(t) is periodic. 
Then, 

N' 

n,n'—0 

= exp{—n{N)) . (23) 

From Eq. ( [T7| ), we know that there exists a constant c 
such that, if n, n' < N' = cN < N, then 


N’ 


n,n'—0 
N' 

- E(d)- 


o-*("+l/2)t 


n—0 


= exp(—n(A^)) . (24) 


Using Eq. (23), it follows that 


N' 




n—0 


= exp(—r2(iV)) . 


(25) 


We now consider the case where = \j). 

We seek to show that {j\U‘^{t)\(f'^) approximates 
(27r/iV)^0(2;^.|{ 7 (i)|^). Equation (22) implies 


N' 


\{j\U^{t)\^^) - = exp(-U(iV)) . 

n—0 

(26) 

Additionally, 


N' 

E' 

n—0 


„-j(n+l/2)t / -I /d 


(Mn) = 


N' 


= (27r/7V)^/^ ^ c„e *(”+^/^^*(a;j|'!/'n) 

n—0 

= (27r/7V)i/4 {xj\U{t)\ip) . 

Then, 

= eM-m)) (28) 


(27) 


in this case as well. 

If |t| > 1, it follows that there exists N = 0(log(l/e)) 
that implies exp(—U(A^)) = 0(e) for both cases of 
or |(^'). Then, N = 0(log(l/e) + N') suffices to pro¬ 
vide an overall precision of order e in the computation of 


Theorem basically shows that scattering amplitudes 
of the continuous-variable QHO can be well approxi¬ 
mated by those of the discrete QHO as long as the di¬ 
mension N scales with the largest quantum number in the 
decomposition of |i^) (and | </?')) in terms of eigenvectors 
of H. Also, N is only logarithmic in 1/e. In particular, 
if the states and can be efficiently prepared 

on a quantum computer, the complexity of computing 
((p''^|C/‘^(f)|(/?‘^) will be dominated by the complexity of 
implementing the evolution operator U^{t) for the corre¬ 
sponding choice of N. 


A. Time evolutions: Trotter-Suzuki approximation 

To compute the propagator on a quantum computer, 
we seek an implementation or simulation of whose 

complexity will be determined by the number of two- 
qubit gates required to approximate the operator. Unfor¬ 
tunately, has large norm, so known results on Hamil¬ 
tonian simulation |151 na [201 EH EH l36] are not very 
useful in the current case; other approaches are neces¬ 
sary. 

Since is diagonal and each entry can be effi¬ 
ciently computed, a quantum computer simulation of 
exp(—i(a:‘^)^t) can be done efficiently. To show this, 
let q{6) = 0(polylog(A/<5)) be the number of two- 
qubit gates required to compute a diagonal entry of 
(x'^)^ within precision <5 > 0. Then, a diagonal en¬ 
try of can be computed within precision e us¬ 

ing q{e/\t\) + 0(polylog(A|t|/e)) two-qubit gates: we 
need 0(log(A|t|/e)) bits (or qubits) to represent (x'^)^ 
within precision e/|t|, and multiplication by t can be 
done using additional 0(polylog(A|t|/e)) gates. Then, 
exp(—i(x'^)^|t|) can be simulated on a quantum com¬ 
puter, within precision e, using 0(polylog(A|<|/e)) two- 
qubit gates. The complexity for the simulation of 
exp(—z(p‘^)^t) is of the same order, since x‘^ and 
are related by the centered Fourier transform and 
exp(—i(p'^)^|t|) = (F/*)“^. exp(—i(x'^)^|t|).F/*. Addition¬ 
ally, , where A is a cyclic per¬ 

mutation [Eq. (|l4|)]. (Note that X^ = II and X^/"^ = 
A“^/^.) An efficient quantum circuit for that re¬ 
quires 0(polylog(A)) two-qubit gates is known p71157] . 

can be implemented on a quantum computer using 
0(polylog(A)) two-qubit gates in a number of different 
ways. For example, X^/"^ can also be decomposed as 
^Af /2 _ p"d_^Ar /2 where Z is the diagonal uni¬ 

tary whose nonzero entries are the roots of unity, i.e., 
exp(i27rfc/A) for k G {0, ...,A — 1}. Then, the diag¬ 
onal entries of are exp(i7rA:) = ±1, resulting in a 

simple and efficient quantum-computer implementation 
of Z^!"^. We also note that when A = 2“ is associated 
with a system of a qubits, X^^^ and Z^^"^ are simply 
the Pauli operators and respectively, acting on 
the first qubit. 

The above results suggest using the Trotter-Suzuki 
product formula to simulate U‘^{t) [TH [T51 [T51 fTOl E5l E5] . 
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Such a formula approximates the evolution operator as 
a sequence or product of shorter time evolutions un¬ 
der and Known results provide an upper 

bound on the number of terms in the formula given by 
M. = 0(exp(6/p)(||i7'^||where 77 > 0 is ar¬ 
bitrarily small, 6 > 0 is a constant, and e > 0 is the 
precision dSl HH- Since ||i^^|| = 0{N), the results 
in [m [ig would imply that the number of two-qubit 
gates required to approximate grows faster than 

being undesirably large in the asymptotic limit. Re¬ 
markably, an improved analysis of the complexity of high- 
order Trotter-Suzuki product formulas allows us to re¬ 
duce the complexity substantially. The basic idea is to 
note that the cost resulting from such formulas actually 
depends on quantities (norms of nested commutators) 
such as (p‘^)^]|| rather than In partic¬ 

ular, a;^, and {a:,p} are a basis of a Lie algebra sp(2) 
of dimension 3. Since our discrete QHO operators ap¬ 
proximate the operators in the continuous-variable case, 
it is possible to show that (p‘^)^, and {x‘^,p‘^} are 

almost a basis of a Lie algebra of dimension 3, when pro¬ 
jected onto the low-energy subspace. This implies, for ex¬ 
ample, ||(5[(a;'^)^, (p‘^)^](5|| = 0{N), where Q is the pro¬ 
jector into the low energy subspace spanned by 
n < cN. In contrast, a simple analysis would have re¬ 
sulted in ||(5[(a;‘^)^, (p'^)^]Q|| = 0{N^). 

The (symmetric) Trotter-Suzuki approximations of the 
evolution operator over a course of evolution time s are 
defined recursively as follows pTlITH] : 

■■= {U^isp)ru^{s - 4sp){U^{sp)r , (29) 

and Ui{s) := Here, p > 

2 is integer and Sp = s/(4 —4^/^^^+^^). For time t, we will 
approximate the evolution operator U‘^{t) by (C/p(s))^, 
where k = t/s and p are chosen to minimize the number 
of exponentials of (a;*^)^ and (p‘^)^ in the product. 

In Appx. B, Lemma we show that there is a 
choice for the dimension of the Hilbert space where 
N = exp{0(-yiog(fV'|t|/e))} -I- 0{N') a choice of p = 
0(A/log(iV'|t|/e)), and |s| = 0(5“^), such that 

||[«(s))^-C/'^(t)]|7A)^)|| = 0(e), (30) 

for all n < N'. We will use Eq. ( [30| to prove the main 
result of this section: 

Theorem 2. Let |(p) = initial 

state of the CV QHO, t = 0(1) the evolution time 
(\t\ ^ o,nd e > 0. Then, there exists N = 

exp(0(Vlog(iV7e))) + 0{N'), p = Q{^\og{N'/e)), and 
|s| = 0(5 P), such that (k = tls) 

m^{s)r-U^{t)]\^^)\\=0{e), (31) 

for all oc ■ The number of exponen¬ 

tials of (x^)"^ and (p“)^ in the product (?7p(s))^ is M. = 
0 (exp( 7 -\/log( 7 V'/e))), where j > 0 is a constant. The 


number of two-qubit gates to simulate {Up{s))^ within 
precision e is M = 0{eyip(ff\J\og[N'/e))), where 7 > 0 
is a constant. 


Proof. The first result is a direct consequence of Eq. ( |30[ ) , 
for the same choices of N, p, and s as in Lemma [^- see 
Eq. (B33l in Appx. B. The only additional error is that 
coming from the fact that the states Itp^) are not exactly 
an orthonormal basis for n < N'. Such an error is order 
vi{N), and the choice of N implies that this error is 
negligible if |t| > 1 . 

The definition of Up implies that the total number of 
exponentials of (a;'^)^ and (p*^)^ in U!^{s) is bounded from 
above by bT [Eq. ([2^]. Also, the choice of s in Lemma 
satisfies 5^1 s| = 0(T). Then, the total number of expo¬ 
nentials in {Up{s))^ is M. = OlfiPk) = 0(|t|5^^), where 
k = t/s. That is, M. = 0{\t\ exp( 7 -y/log(fV'|t|/e))), for 
some constant 7 > 0 . 

To achieve overall precision e, each exponential of 
(x'^)^ or (p'^)^ in (Up(s))^ has to be simulated within 
precision Oie/M). Then, the above discussion on the 
cost of implementing diagonal unitaries implies that 
the number of two-qubit gates for each exponential is 
0(polylog(A^|s|AI/e)) or 0(polylog(A^|t|5^/e)). [Note 
that [spl = 0(s) and (s — 4sp) = 0(s).] The choice 
of N and p imply that log(Af|t|5P/e) = 0(log(|t|/e)) -I- 


0{^/\og{N'\t\/e)) -\- 0(logN'), and we can safely bound 
this quantity by 0(log(Af'|t|/e)). Then, the total num¬ 
ber of two-qubit gates is 0{M polylog(A^'|t|/e)). That 
is, there exists a constant 7 > 0 such that the number of 
two-qubit gates to simulate {Up{s))^ within precision e 
is AI = 0{\t\ exp( 7 -y/log(A''|f|/e))polylog(W|t|/e)). We 
let X > 1 and fc > 0 a constant. Then, there exists a 
constant 7 such that < e^^. It follows that 

TW = 0{\t\ exp( 7 -yiog(A'|t|/e))) and the result is ob¬ 
tained setting t = 0 ( 1 ). □ 


We illustrate the results of this section with several nu¬ 
merical simulations. In Fig. [^we show the error depen¬ 
dence of the Trotter-Suzuki approximation as a function 
of n for fixed s and p. A worst-case analysis indicates 
that this error would be 0 (n^^^’+^)/^) [TKl 1 ^ . However, 
a linear scaling in n is observed, for n < N/2. This is 
a main reason for the quantum speedup; see Appx B. In 
Fig. [5] we choose p and s according to Thm. [^and show 
that the approximation error is indeed bounded by e. 


B. Subexponential-time quantum algorithm for 
computing scattering amplitudes 

We can combine Thms. and to construct a quan¬ 
tum algorithm to simulate the QHO and obtain the de¬ 
sired propagators {ip'\U{t)\(p). The quantum algorithm 
has three basic steps: i- the preparation of an initial 
state, ii- the implementation or simulation of the evo¬ 
lution operator, and iii- a projective measurement. The 
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FIG. 4: The difference between U‘^{s)\il>n) and Up{s)\'ipn) as 
a function of n, for N = 400, p = 4, and s = 1. The scaling 
is almost linear in n, as suggested by Lemma in Appx. B. 
Additional simulations also observe a linear scaling in n for 
higher values of p. 



FIG. 5: The difference between and {Up{s))’°\'ipn) 

as a function of the dimension, for t = 7r/2 and n = N/2. 
The values of p and s where set as determined by Thm. 
That is, according to Lemma in Appx. B, we chose p = 
rVlog((n + 2)t/e)/(21og(5))l, k = \t/{e/{{n + 
and s = t/k, for e = 0.001. Note that with these choices, 
s = exp(—0(\/log(Aft/e))) and the number of terms in the 
product (t/p(s))*^ is A4 = = 0{t exp( 7 y^log(Nt/e))), for 

some 7 > 0. The error induced by the Trotter-Suzuki ap¬ 
proximation in this case is smaller than e. The “jump” at 
N « 500 corresponds to an increase of p from 2 to 3, due to 
the increase in N. 


desired expectation value can be obtained within arbi¬ 
trary accuracy after repeated executions of these steps. 
The simulation of the ev olution operator induced by 
was discussed in Sec. 


IV A 


We let and W^, 
be the N x N unitary matrices that prepare the states 
and from e.g. |0), respectively. We also let 

so the propagator of interest is 
(0| V‘^(t)|0). To measure such expectation values at pre¬ 
cision e, we can implement the quantum circuit in Fig.[^ 
0(1/^/e) times [^lEO]. Otherwise, we can implement the 
quantum methods described in [31] for optimal quantum 
measurements of overlaps, where the number of repeti¬ 
tions can be improved to 0(l/e). 

In the rest of this section we assume that the preci¬ 
sion e <C 1 is constant, i.e., e = 0(1). We also assume 



(0|F’‘(t)|0) = V+) 



FIG. 6: Quantum circuit to compute the propagator 

(0| V‘^(t)|0), where V‘^{t) is unitary [S] [40]. The filled cir¬ 
cle denotes the controlled operation on the state |1) of the 
ancilla qubit; that is, the corresponding unitary operation is 
1 <8) |0)(0| + (8) H denotes the Hadamard gate 

that maps the state of a qubit as H|0) = (|0) -I- \ l))/\/2 
and H|l) = (|0) — \l))/\/2. The single qubit operator cr'*’ 
is (Tj; +i<Jy, with a a the Pauli operators. Since cr"'" is not Her- 
mitian, the computation of its expectation value can be done 
by repeated projective measurements of and Oy indepen¬ 
dently (i.e., with repeated executions of the circuit). 


that and W^, can be efficiently implemented using 
0(polylog(V)) gates. Then, Thms. and [pimply 

Theorem 3. Let \ip) = Y.n=o<^n\ipn), W) = 

'Ylin=o^'n\i’n) or \ip') = (27r/iV)^/^|xj) and t the evolu¬ 
tion time. Assume |t| > 1 and, with no loss of gen¬ 
erality, t = 0(1). Then, there exists a quantum algo¬ 
rithm Q that outputs {ip'\U(t)\ip) at arbitrary accuracy. 
Q requires M.q = 0 (exp( 7 '-\/log(W))) two-qubit gates, 
where > 0 is a constant. 

Proof. Let N = exp(0(y/log(iV'|t|))) -I- 0{N'), p = 
0{\/\og(N^), and |s| = exp(-0(v/log(fV'|t|))) as in 
Lemma or the proof of Thm. for e = 0(1). The 
subadditivity property of errors implies (fc = t/s) 

{p'\U{tM = {p'^\{U^{s)r\p<^)+0{e) 

= (0|(W^^,)^(C^p W)'=<|0) + 0(e) , (32) 

with e ^ 1. Then, we can use the circuit of Fig. 
a constant number of times, with a unitary V'^{t) = 
{Up{s))'^W^ to output the desired propagator. 
Since W^, and can be efficiently implemented, the 
complexity of the quantum algorithm is dominated by the 
number of two-qubit gates needed to implement (C/p (s))^ 
within the desired accuracy. Theorem [^implies that such 
a number is 0 (|t| exp( 7 'A/log(iV'|t|))), for some constant 
7 ' > 0. The result follows by assuming, with no loss of 
generality, t = 0 ( 1 ). □ 

It is important to remark that the complexity of the 
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quantum algorithm satisfies 

lim Mq/N'^ = 0, (33) 

N'—^ao 

for all ry > 0, and lim^v'-j-oo log(-/V')/A^Q = 0. That is, 
the complexity of the quantum algorithm is subexponen¬ 
tial in log{N'). Since classical algorithms are expected 
to require poly(-/V') operations to compute the propaga¬ 
tor in the worst case, our quantum algorithm provides a 
superpolynomial quantum speedup. 


Proof. First, we choose 6 = 0(log(l/e)) so that the right 
hand side of Eq. (34) is 0(e). In the large N limit, |'0o) 
can be safely replaced by |(/)q) in Eq. (34), as the error 
of this replacement is exponentially small in N and thus 
negligible. Next, we note that we can approximate 
within error 0(e) by 


jo 

(X ^ \j) , (37) 

j = -jo 


V. EIGENSTATE PREPARATION 

We now investigate ways of simulating and prepar¬ 
ing low-energy eigenstates of via the action of uni¬ 
tary operations acting on simple initial states. In part, 
this section addresses the second goal of a QS, namely 
the computation of expectation values on various eigen¬ 
states of the QHO, which can be obtained using the 
techniques presented in previous sections if we replace 
the initial state |(/?) by the corresponding \ijjn) (or 
by \(j)n))- The results of this section may be of inde¬ 
pendent interest; e.g., quantum algorithms to prepare 
states with Gaussian-like amplitudes are important in 
other cases mi- 

We first focus on the premration of the ground state 
|(/)q). In Appx. C, Lemma we prove in the large-A 
limit, 

lll^d^ _ e—= 0(exp(-fl(<5))) , 

(34) 


where the initial state is 

N/2 

(1/v^) exp(-jV(2(5)) \j) , (35) 

j=-N/2 

and (5 > 0. The constant k is for normalization pur- 
porses and a(t) is an irrelevant global phase that can 
be computed exactly. The evolution times satisfy t = 
•\/(7^(2 — 4(t 2)/2 and t' = l/(4t-|-4(T^/t), and = ttS/N. 
This result was obtained by realizing that in CVs, the 
quantum state I'i/^o) can be obtained from an initial state 
with Gaussian-like amplitudes by evolving with the free- 
particle Hamiltonian (i.e., —p^). The result follows by 
approximating the GV case after a proper discretization. 

Lemma |9] allows us to state the first result of this sec¬ 
tion. 

Theorem 4. Let e > 0. Then, there is a unitary W‘^ 
that satisfies 


10)11= 0(e) (36) 

in the large N limit. can be implemented on a quan¬ 
tum computer using a number of two-qubit gates that is 
polynomial in log{N/e). 


with jo = 0(-\/log(l/e)<5) = 0(log(l/e)). The state of 
Eq. (37) can be prepared with complexity polynomial 
in log(l/e) using standard techniques. We write for 
the unitary that prepares such a state and define = 

the choices of t and t' given 
above. Because t = 0(1/VN) and t' = 0{\fN) in the 
large N limit, the unitaries ^ * and e*0 ) * can be 
implemented on a quantum computer with complexity 
polynomial in log(l /e) an d log(A) [i.e., polynomial in 
log(iV/e)] - see Sec. ra □ 


In Fig. [^we show the exponential decay of the error as 
a function of 5, as stated by the theorem. We note that 
quantum methods to prepare states with Gaussian-like 
amplitudes were also proposed in mnn]. 


0.0020 


— 0.0015 

0.0010 

■a o 

0.0005 


1 2 3 4 5 5 

FIG. 7: The norm of the difference between the state |jio) 
and the evolved state |yi‘^(t)) = ) t 

as a function of S and for N = 800. Numerical simulations do 
not show significant changes for larger dimensions. 


To prepare the other eigenstates |()>)[), with n> 1, we 
define a discrete version of the Jaynes-Gummings (JG) 
model: 


Hjc = G CTa; - (g) , (38) 

where Ua are the corresponding Pauli operators acting 
on a the Hilbert space of an ancillary qubit. In GV, the 
evolution induced by the JC model eventually transforms 
the state lipn) |0) into |'0ra-i-i) |1)) providing a unitary op¬ 
eration to prepare eigenstates of the QHO from Itpo)- We 
will show that something similar occurs in the discrete 
case. 
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For n > 0, we define the normalized states 

K±) = ^[|-/>^)|0)±|</>^+i)|l)]. (39) 

These are approximations of the eigenstates of Hf^. In 
Appx. D, Lemma we show that if n < < cN, for 

some constant c > 0, 


\\{HfcTV^)hi±)\\=MN). (40) 

We use Eq. p0| to prove the second result of this section. 

Theorem 5. Let e > 0 and tn = TTl{2y/n + 1). Then, 
there exists N = 0(log(l/e) + N') such that 

||e-^.V.|0d)|o)+,|^d^^)|l)||=O(,)^ (41) 


for all n < N'. 


Proof. We let N > N'/c. Then, Eq. (40) implies 


_ p-i(±Vn+ltr 


H\li±)\\ = \tn\MN), (42) 


for all n < TV'. For = t :+ 1), this implies 
±i)|^d^)|| = jyi(iV), and then 


|0) |l) II = ^,(iV) . (43) 


Then, there is = 0(log(l/e)) so that the right hand 
side of Eq. (43) is 0(e}. □ 


We can combine Thms. and to prepare approxi¬ 
mations of other eigenstates |(^((), n > 1, by a sequential 
action of (I® on the initial state |())q) |0). 

The Pauli operator ctj, is necessary to transform |1) —>■ |0) 
for the state of the ancilla qubit at each step. Then, 


n—1 

lO |0) « n ® |<(>^) |0) . (44) 

n'—O 


We now seek a quantum algorithm to prepare the 
eigenstates \4>n), up to a given approximation error. This 
requires giving a quantum circuit to approximate each 
in Eq. (44|). Since the unitaries 


and <Si<yv)t ^^n be simulated within precision e us¬ 

ing a numb er of two-qubit gates that is polylog(A^|t|/e) 
(Sec. IVA), we will use the Trotter-Suzuki approxima¬ 
tion. In Appx. D, Lemma EH we show that if s = 
0[el\Jn P 1) for some n < N' < cN, then 

||jg-i(a;‘‘®(T,„)s/V2gi(p‘‘®(T„)s/V2 _ ^|)|| = 0{e^) . 

(45) 


The proof uses a simple Trotter-Suzuki approximation 
and the scaling with e can be improved using higher order 
approximations. We can use this result to prove: 


Theorem 6. Let e > 0. Then, there is a quantum circuit 
that satisfies 

l|VE"|^^)|0)-|</.)()|0)||=O(e) (46) 

for any given n < N' < cN, where N is the dimension 
of the Hilbert space and c > 0 is a constant. can 

be implemented using a number of two-qubit gates that is 
0((n2/e)polylog(V/e)). 

Proof. The quantum circuit is 

n—1 

n'—O 

(47) 




Here, m = 0{n/e) so that each term in IF‘^ introduces 
an error Oiejn') as implied by Lemma 11 Because we 


work in the asymptotic limit, approximation errors that 
are exponentially small in N or \/]V are negligible. Then, 


|0) - |0) II = 0{{e/n)n) = 0(e) . (48) 


The number of terms in the product is 0{n'^le). The 
number of two-qubit gates is then 

0((n^/e)polylog(V/e)) . (49) 


□ 


VI. THE QUARTIC POTENTIAL 

We now analyze quantum algorithms to simulate the 
evolution operator of a quantum system with Hamilto¬ 
nian H — 4(p^ -(-x"^). We will use the same discretization 
as that for the DQHO, where and where defined in 
Sec. |HI| Then, 

H^ = li{p‘^r + {x^r). ( 50 ) 

In contrast to previous sections, we only conduct a nu¬ 
merical analysis here and state some observations from 
numerical results. In part, this is due to not having an 
exact solution in this case. Our simulations suggest a 
polynomial speedup for the computation of scattering 
amplitudes. 

In Fig. we plot the eigenvalues E!^{N) of 
for different dimensions N and as a function of n = 
0,1,..., V — 1. Taking dimension 4000 as a reference, 
in Table U we look for the maximum value of n such that 
I ^^(V) —111)^(4000) I is below some threshold e* <C 1. The 
values of the ratios r* suggest that the eigenvalues of a 
large sector of the low-energy subspace of the CV system 
can be well approximated by those of the discrete sys¬ 
tem. While P does seem to decrease in N, the scaling 
does not seem to be of the form 1/N^, for some y > 0, 
but rather of the form 1/log V. If this is the case, then. 
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to approximate up to the n-th eigenvalue of H, it suf¬ 
fices to choose N = 0{n\ogn) for the discrete system. 
Nevertheless, both the classical and quantum algorithms 
to simulate this system will have complexity that depend 
on the same value of N, and the dependence of r® on N 
is unimportant to demonstrate a quantum speedup. 



FIG. 8: Eigenvalues of H'^ for N = 4000 (blue), N = 3200 
(yellow), and N = 1600 (purple). 


N 





100 

16 

0.16 

12 

0.12 

200 

30 

0.15 

27 

0.135 

400 

57 

0.1425 

54 

0.135 

800 

105 

0.1312 

99 

0.125 

1600 

181 

0.1131 

177 

0.1106 

2000 

216 

0.108 

212 

0.106 

3200 

314 

0.0981 

310 

0.0968 


TABLE I: n', maximum value of n such that |i?)5(N) — 
£’(((4000)1 < for = 10“® and ^ = 10~^. The ratios 
are r* = n*/N. 


We can use the Trotter-Suzuki approximation to sim¬ 
ulate the evolution operator = exp{—iH^t}. This 

approximation splits the evolution operator into a prod¬ 
uct of exponentials or short time evo lution s under (p'^)^ 
and (x'^)^ and, using the result of Sec. IV A[ each of these 
exponentials can be simulated efficiently. In contrast to 
the QHO, the operators x'^ and do not form a small di¬ 
mensional Lie algebra and tight error bounds from high- 
order Trotter-Suzuki approximations may be difficult to 
obtain. The recursive definition for Up{s) is given in 


Eq. (29), assuming 


Utis) := e- 


(51) 


A worst-case analysis of the Trotter-Suzuki formula 
results in an approximation error bounded by ep(s) = 
0(|s|2p+i||iLd||2p-ri) = 0((|s|V2)2p+i) [m [Ill [IHl [m 
1551 155] . However, as in the case of the DQHO, we 
would expect that the error for the current case is sig¬ 
nificantly smaller than that for the worst case. This is 
because the operators in H, while they do not form a 
finite dimensional Lie algebra, posses an algebraic struc¬ 
ture that results in an effective norm for nested com¬ 


mutators that is significantly smaller than the prod¬ 
uct of the effective norms [32]. In Fig. we plot the 
error ||(17p(s) — tf'^(s))|^())()|| as a function of n, for 
|s| = 0(1), and p = 1,2,3. Here, |(/)^) are the eigen¬ 
states of in Eq. (50) and cannot be approximated 
by the discrete Hermite states \tpn). The results sug¬ 
gest \\{U^{s) - U‘^{s)) lO II = 0(|s|2p+in(2p+4)/3). The 


order dependence in s follows from the analysis of the 
high-order Trotter-Suzuki approximations and was veri¬ 
fied by additional numerical simulations, to assure that 
we are in a region of convergence. 




FIG. 9: The error from high-order Trotter-Suzuki approx¬ 
imations (red dots) as a function of n. The results are for 
dimension N = 1600, s = 0.01, and for n = 0,1,..., 181, 
according to Table jl] Our simulations suggest that the er¬ 
ror is 0(s2p+^nP ), where p' depends on p. (a) p = 1. The 
blue line indicates a fit with the function /i = 3.44.10”^ . 
(b) p = 2. The blue line indicates a fit with the function 
/i = 7 . 8 . 10 “^ 2 j 28 / 3 , (c) p = 3. The blue line indicates a fit 
with the function /a = 3.3.10“^^ 


Under these numerical observations, we can analyze 
the complexity of a quantum algorithm that computes 
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scattering amplitudes 


{ip'\ U{t) \ip) , 


(52) 


within some precision e. For simplicity, we assume 
that the initial and final states, \(p) and |(/?'), can be 
written as linear superpositions of the eigenstates \4>n), 
with n < N'. This suggests that we can approximate 
Eq. (521 by , where the dimension is 

N = 0{N' \og{N')), and and are quantum 

states obtained by replacing |^„) by |(/)^) in the spectral 
decompositions of \^p) and |(/?'), respectively. We assume 
that and can be efficiently prepared so that the 
main cost of the algorithm is that from simulating U^{t). 
We split the evolution time into k parts of size s = t/k. 
Using the {2p + l)-th order Trotter-Suzuki approxima¬ 
tion, our numerical simulations suggest that the error is 
bounded by kuj{ui' |s|)2p+17V(^p+^)/^. uj and w' > 1 are 
constants. Replacing s by t/k, the number of terms in 
the product formula is 


depend on N' and the precision parameter e. In partic¬ 
ular, for constant precision, we will assume that N = 
0(poly(fV')), an assumption that is satisfied by a large 
class of one-dimensional quantum systems such as the 
QHO or the quartic potential. 

Theorem 7. Let e > 0 and t be the evolution time, |t| = 
U(l). Then, there is a quantum circuit that satisfies 

\\{W‘^-U‘^{t))\^^)\\<e, (56) 

for all cx Here, Ifi//} are the eigen¬ 

states of and U‘^(t) = exp{— can be im¬ 
plemented using 0(A^|t|polylog(A^|t|/e)) two-qubit gates. 

Proof. We use our results in |21j for Hamiltonian simu¬ 
lation and do some modifications to obtain the desired 
result. In |2I] we showed that for a d-sparse time depen¬ 
dent Hamiltonian A{t) acting on q qubits, the evolution 
operator can be approximated within precision e using 


k5P = 0 


/|^|l-Hl/2p^(l-H2/p)/3 

i ei/2p 



(53) 


0(gd^|jH||max|i| 


l0g((d^PI|max+||i||)|t|/e) 

log log(d2||H||j„ax|i|/e) 


(57) 


As p grows large, the number of terms can be made 
for arbitrary small r]. Then, the 
quantum circuit that approximates the evolution oper¬ 
ator U^{t), needed to compute Eq. (52), can be imple¬ 
mented using a number of two-qubit gates that is 


' |f|l+P/Vl/3+4p/3 

Mq = 0{ - - - -polylog(Aft/e) 


(54) 


This complexity represents a polynomial quantum 
speedup, with respect to N, over the quantum algorithm 
that approximates scattering amplitudes for the quartic 
potential. Classical algorithms, in the worst-case, may 
require computing and obtaining the spectral properties 
of U^{t), which can be done in complexity 0{N'^), for 
cr > 2. 


VII. ONE-DIMENSIONAL QUANTUM 
SYSTEMS: AN UPPER BOUND 

We now present an upper bound on the complexity 
of simulating the evolution operator of one-dimensional 
quantum systems described by 77 = -\- V{x), where 

V(x) is the potential, i.e., some operator that depends on 
X. If \4>n), n = 0,1,..., denote the eigenstates of H and 
\fi) is the initial state, we will assume \ fi) = \4‘n) 

and ||U(a;) \4in) || = 0(poly(7V')), for all n < N'. 

We also assume that there exists N < oo, the dimen¬ 
sion of the Hilbert space, such that, if 

H'^ = \ip‘^r + V{x<^) , ( 55 ) 


two-qubit gates, where A = dtA and the norms are the 
maximum norms in the time interval [0,t]. The gate 
cost results from a decomposition of Aft) in terms of 
0(d^||A||niax) unitary operators. Since is not sparse 
and its norm can be larger than 0{N), Eq. ( |57| ) would 
result in a large gate complexity in this case. To overcome 
this difficulty, we analyze the evolution operator in the 
interaction picture. We then define the time-dependent 
interaction Hamiltonian 

Hf{s) = leiV{xfis^pd-^2^-iVixfis 

= , (58) 


and denote Uf{t) for the corresponding evolution. Since 
is 1-sparse and ||a:'^|j = 0{'/N), we obtain ||77f ||max = 
0{N) and 77/(s) can be decomposed as a sum of 0{N) 
unitary operators. In addition, ||77/|| = 0(poly(7V)) un¬ 
der the assumptions and q = 0(log N). This implies that 


the gate complexity to simulate Uf as given by Eq. (|57|) 


is 


0{N\t\log^{N\t\/e))=d{N\t\) 


(59) 


if |t| = Q(l). We note that U<^{t) = and 

that can be efficiently simulated with complex¬ 
ity polylog(V|7|/e) as explained in Sec. IV A □ 


VIII. THE FRACTIONAL FOURIER 
TRANSFORM AND RELATED WORK 


the scattering amplitudes of the CV system can be well The evolution induced by the QHO results in a trans- 
approximated by those of the discrete system. TV will formation referred to as the fractional Fourier transform 
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(frFT), which corresponds to an arbitrary rotation in 
phase space. The frFT has been proven useful in sig¬ 
nal analysis |43H45j . in noise filtering in particular, when 
the noise does not have a well defined frequency spec¬ 
trum. The evolution operator U^{t) can then be inter¬ 
preted as an approximate version to a “discrete” frFT, 
and our results may prove useful in the design of classical 
or quantum algorithms for discrete signal analysis. This 
would require efficient methods for encoding and decod¬ 
ing of signals, which may exist under some assumptions 
such as sparsity. 

The split-step (Fourier) method is a well known tech¬ 
nique to solve the nonlinear Schrodinger equation in 
quantum mechanics and in fiber optics (c.f., |46jb As 
in our case, the idea is to evolve the initial state ac¬ 
cording to small step evolutions under the correspond¬ 
ing operators. This also requires a discretization of the 
continuous-variable coordinates. Our results on Trotter- 
Suzuki approximations can then be used and generalized 
to bound the errors induced by the split-step method. 

The development of quantum simulation methods for 
continuous-variable quantum systems, including quan¬ 
tum chemistry, is an active area of research (c.f., [71151 ITT] 
and references therein). Commonly, the complexity of 
such methods is polynomial in the energy of the system. 
As an example, in |29j , a quantum algorithm for approx¬ 
imating the ground state energy of a continuous-variable 
quantum system is provided. The algorithm works under 
some assumptions on the energy potential and its com¬ 
plexity is proportional to d, the number of state variables 
(or particles). Another example is [30], which presents a 
quantum algorithm to compute scattering probabilities 
in a certain quantum field theory theory). The com¬ 
plexity of such algorithm is also polynomial in the energy. 
In contrast, our quantum algorithm to simulate the time 
evolution operator of the QHO has complexity that is 
subexponential in logA^', where N' denotes the relevant 
energy scale of the problem. The QHO provides a ba¬ 
sis for the quantization of the electromagnetic field and 
quantum field theories, and our algorithms are expected 
to find wide applications in quantum simulation. 

IX. CONCLUSIONS 

We provided a quantum algorithm to approximate the 
propagator of the QHO within arbitrary accuracy. For 
precision e > 0, the complexity of the algorithm is A4q = 
0(exp(7-\/log(A^/e))), where the evolution time can be 
assumed to be constant, and N is the relevant energy 
scale of the simulation. Asymptotically, Mq/{N)^ —>• 0, 
for any 77 > 0, so the complexity of the algorithm is 
subexponential in log(Af). Remarkably, this represents 
a superpolynomial speedup over the corresponding clas¬ 
sical algorithm to compute the propagator, whose com¬ 


plexity is 0{N). Our results consider a refined analysis 
of the error of high-order Trotter-Suzuki approximations. 
This analysis works in this case because the operators 
under consideration form a Lie algebra of dimension 3 
[i.e., sp(2)]. We can then use properties of commutators 
to show that the error induced by a high-order Trotter- 
Suzuki approximation is significantly smaller than that 
obtained considering the worst-case scenario; recent re¬ 
sults consider this problem more generally [32] . Our 
quantum algorithm considers a discrete version of the 
QHO whose low-energy spectrum can be shown to repro¬ 
duce the properties of the QHO with very high accuracy 
(i.e., the approximation errors decay exponentially with 
N). 

We also provided quantum algorithms to compute 
spectral properties by preparing approximations of the 
eigenstates of the QHO. These algorithms have complex¬ 
ity polynomial in log(A^)/e and may be of independent 
interest (e.g., the ground state has Gaussian-like ampli¬ 
tudes). Here, N is the number of points in the discretiza¬ 
tion or dimension of the Hilbert space. To prepare such 
states, our quantum algorithms simulate the evolution in¬ 
duced by a version of the Jaynes-Cummings model that 
is used in quantum optics. This evolution is also approxi¬ 
mated using a high-order Trotter-Suzuki approximation. 

Last, we presented quantum algorithms to simulate 
more complex one-dimensional quantum systems. For 
the case of a quartic potential, we presented numerical 
evidence for the existence of a method that simulates the 
evolution operator with complexity Q = 0(Ai^/^+°^^^), 
if t and e are constant. This method is also based on high- 
order Trotter-Suzuki approximations. Our result repre¬ 
sents a polynomial quantum speedup over the classical 
method in this case. The quantum advantage may be 
a result of the algebraic structure satisfied by the op¬ 
erators in the Hamiltonian [32]. We also showed how 
general quantum systems can be tackled on the basis of 
recent results in |2I] that simulate the evolution opera¬ 
tor by implementing its Taylor series decomposition. We 
proved an 0{N) bound for the complexity of simulating 
the evolution operator under fairly general assumptions. 

We conjecture that some of our results can be general¬ 
ized to provide subexponential time quantum algorithms 
to simulate the evolution of other continuous-variable 
quantum systems. 
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Appendix A: Properties of the discrete QHO 

Most of our results can be obtained from approxima¬ 
tions of integrals appearing in the continuous-variable 
case as finite sums appearing in the discrete variable case. 
For completeness, the Hermite functions are 

ip„{x) = ^ , (Al) 

^/2^n\y/TT 

where Hn{x) is the (physicists’) n-th Hermite polyno¬ 
mial, n > 0, 

Hnix) = (-l)"e- . (A2) 

The orthogonality (or normalization) condition is 
/dx = Sn,m, and a useful property is 

xilJriix) = + l/(n-b l)/2V'n+l(x). In 

Dirac’s bra-ket notation, these are {tpm\4’n) = dn,m and 
xji/'n) = \/ri72|V’n-i) + \/{n + 1)/2|V'„-hi), respectively. 
The Hermite functions are eigenfunctions of the Fourier 
transform with eigenvalues (—*)", n > 0. Then, we can 
write F for an operator that applies the Fourier trans¬ 
form and Fjf/’n) = (-i)”|'0n). 


Unless a region of integration is given explicitly, we 
assume x S (— 00 , 00 ). Similarly, sums are infinite unless 
otherwise stated. In all cases, I > 0, k, and n > 0 are 
integer numbers. N > A denotes the dimension of the 
Hilbert space for the discrete QHO. 

We will show that our main results follow from state¬ 
ments about the “tails” of the Hermite functions for 
X ^ [—kT/2,kT/2], where T = y/2'KN and fc > 1. For 
even n, the Hermite polynomials satisfy 

(A3) 

and there is a similar upper bound for odd n (see in¬ 
equality 8.954 of [IH]). Stirling’s approximation implies 
ci\/ri(n!e)'^ < n! < C 2 y/n{nje)'^, with 0 < ci = < 

C 2 = e. Then 

IV’n(x)| (A4) 


for some constant C 3 « 0.7 and x > 0. In the event that 
x/2 — > P'/N, for some P > D, Eq. (A4| implies 

ipn{x) = 0{exp{—xpVN)). A similar result is obtained 
for odd n. Also ([45]), 


|iF„(x)| < C 4 V^ 2 ”Aea: /2 ^ 


(A5) 


which implies 


IV’n(x)| < 1 , (A 6 ) 


where C 4 « 1.0864. 

The results described in Secs. IIII Al and IIVI will fol¬ 
low from the following lemmas, which suffice to provide 
analytical proofs. For simplicity, we use vi{N) to de¬ 
note the order of a function that decays exponentially 
in N. That is, vi{N) = exp(—H(A)), and there ex¬ 
ists a constant /? > 0 such that, if f{N) = vi{N), then 
f{N) < exp(—/3A^). Note that (A^“ exp(—H(A1))) and 
i/i(A^)“ are also exp(—H(A)) for any constant a > 0 
and sufficiently large N. Then, if f{N) = N°‘vi{N) or 
f{N) = we have f{N) = vi{N). The value of 

the constants appearing in these lemmas, such as p, can 
be improved with more detailed analyses. 


Lemma 1. Given I > 0, there exists a constant c > 0 
such that, for all n < cN and all k > 1, 


dx (V’n(x))^x' 


lkT/2 




(A7) 


Proof. We will choose c > 0 so that x/2 — > T/4 — 

> P'/N, for some ,3 > 0 , in the integration region. 
For example, we can choose c < 7r/16 and P to satisfy 
P = {y/TTj& — y/2c) > 0. Then, Eq. (A4) implies ifnix) = 
0{eyi.-p{—xP'/N)) for x > T/2. If is the integral of 
interest, explicit calculation gives 


h = 0{exp{-kTpVN){kVNy) . (A8) 

Since I = 0(1), T = H(ViV), and (ky/NY) = 
exp(Z \og{ky/N)), we obtain = (i/i(A^))^. 

□ 
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We note that the value of jS in the exponential of 
Eq. (A 8 ) depends on c (or n). For example, if we 
are interested in the case where n = 0 , we can choose 
P = Y^ 7 r /8 so that Ik = 0(exp(—fc(7r/2)iV)(fc-\/]V)0- 


Lemma 2. Given I > 0, there exists a constant c > 0 
such that, for all n,m < cN, 


- {f^mW\'tpn)\ = Vl{N) . (A9) 


Proof. From the property of the Hermite functions, we 
note 


i£y\^n)= E ^ = E 

n'<n+Z n'<n+Z 

(AlO) 


where |c„'| = Since I = 0(1), the lemma follows 

from showing that (V'mIV'n) approximates Sn^m within 
precision that is exponentially small in N. 

The remainder of the proof has two parts. First, we 
will show that 

{27r/Nf^Yl )'^n 

3 


with Xj = jy^^TiTiV, approximates (^^l^n) 
sired order. The Cauchy-Schwarz inequality implies 


oo 


I E < 

j=Nl2 


- E E '^niXj) 


\i=N/2 


\j=N/2 


(A12) 


In this case, Xj > \pKWj2. We will choose c so that 
Xjl2 — >T/4, — \/2n > jSyfN, for some /9 > 0. That 

is, c and /3 can be those of Lemma[^ and n < cN. Then, 
Eq. (A4) implies tpnixj) = 0(exp(—XjjdvW)) for xj > 
T12 or, equivalently, for j > N/2. Explicit calculation of 
the sum implies 


In the second part of the proof, we will show 
that Eq. ( [All ) also approximates the expectation 
ii’mli’n) = (>n,rn- We use the identity ^(x - Xj) = 
(A^/(27r))^ /^ g-jfcTa;^ where S(x) is the Dirac delta, to 
write Eq. ( |All| as 

j dx V'm(a;)V'n(a;)e“*''^“ = ' 

k k 

(A15) 

The term with /c = 0 is {'ipmlf’n), so we need to show that 
the sum of the terms with fc 7 ^ 0 is small and satisfies the 
desired bound. The Hermite functions are also eigen¬ 
states of the Fourier transform and ('i/'m|e“*^^'^ \ipn} = 
(—|'0n), where we used Fx{F)^ = p. 
Note that is the space translation operator. Then, 

each term of Eq. (A15) can be written as 


dx 'ljjrn{x)'lpn{x + UT) 


(A16) 


because exp{—kTdjf}'fin{x) = ^/’„(x -|- kT). 

We are then interested in showing that 
/ dx 'iprn(x)'ifn{x F kT) is Small when k 0. From 
symmetry arguments, it suffices to analyze the case 
k > 1 only. We write i>m{x) = 'ifm{x) + (t>m{x), 
where ipmix) = ipm{x) if x < kT/2 and ipmix) = 0 
otherwise. We use a similar decomposition for 
ifnix) = i’nix) Pnix), where ipn{x) = tpnix) if 
X > kT/2 and 'ifn{x) = 0 otherwise. It follows 
that ifn{x)iprn{x) = 4>m{x)4>'nix) = 0 because these 
functions are supported in disjoint regions. Also, 


Eqs. (A4| and ( |A 6 | ) imply 'fn{x)4>m(x) = and 

(p^{x)tpm{x) = vi{N)^ so that Eq. ( |A16| ) is also of order 
vi{N)^. Summing over fc 7 ^ I implies that Eq. (All) 
can be approximated by (' 0 m|' 0 n) = ^n,mj the term with 
fc = 0, or by (t/’mlV'n)) within precision J^i(iV). □ 


The previous analysis implies: 

Corollary 1. There exists a constant c > 0 such that, 
for all n,m < cN, the discrete Hermite states are almost 
orthonormal: 


{2t:/N)^/^^ ■ifl{xj) = vi{N) , (A13) 

j=N/2 


for sufficiently large N. This coincides with the result of 
Lemma H] when fc = 1. The same result can be obtained 
if we replace n —>■ m, with the assumption that m < cN. 
Then, from the Cauchy-Schwarz inequality, 


^ E ^rniXj)'llJn{Xj) 


j=N/2 


= MN) , 


(A14) 


impl ying that (V'mIV'n) can be approximated by 
Eq. (All) within precision that is exponentially small 
in N. 


l(V'mlV^n) -<5n.m| = • (A17) 

Proof. It is a direct consequence of Lemma for I = 

0. □ 


For all n' > 0 integer, we define the states 



with 

00 

fpn'iXj) = E ’^ri'iXj + kT) . 

k— — oci 


(AI 8 ) 


(AI9) 
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Such states will be useful to prove the following lemmas. 
Remarkably, the are eigenvectors of the centered 

Fourier transform with eigenvalues (—*)” ■ To show this, 
we work in the bra-ket notation and let \Cj) = \xj + 
kT) be the states that represent the corresponding Dirac 
combs (sums of Dirac deltas): {Cj\il)n')='4’n'[xj). Then, 






n) 


XI b') (Cj#"') ■ 

3=-N/2 


(A20) 


The properties of the Fourier transform when acting on 
the Dirac comb implies 

N/2-1 

F\C,) = {l/y/N) X . (A21) 

j'=-7V/2 

The centered Fourier transform has a similar action on 

|j): 

Af/2-l 

= (l/TiV) X . (A22) 

3'=-N/2 


Then, 


^ ^ 3=-NI2 

(A23) 


where we used F\tpn') = (—0" l^n')- 
We note that 

Af/2-l 

lll^n) - l^n)f = (27r/^)^^^ X \i’n{x 3 ) - ^n{Xj)\'^ 

3=-NI2 

N/2-1 

= {2TT/Nf/^ X \Y.^n{x,+kT)\'^ . 


3=-N/2 fc/0 


For A: 7 ^ 0, we obtain \xj + kT\ > T/2. Then, if c and /3 
are the constants in Lemmaj^and Lemmaj^ and n < cN, 
Eq. (A4) implies \ijjnixj + kT)\ = 0{exp{—xP\/N)), with 
X = \xj + kT\. Consider, for example, fc > 1. In that 
case, \xj + kT\ > kT/2 and explicit calculation gives 
X]fc>i IV'n(a;j + kT)\ — 0{vx{N)) for sufficiently large N. 
The result for fc < — 1 is similar. It follows that 




(A25) 


Lemma 3. Given I > 0, there exists a constant c > 0 
such that, for all n,m < cN, 

|(V'ml(/)'lV'n) - (V'mlp'lV'n)l = ■ (A26) 

Proof. Because approximates I'tpff) for n < cN, 
where c is the constant of Lemma we can transform 


the operators p and with the corresponding Fourier 
transformations, and use Lemma to obtain the desired 
result. Since lip'll] =_||x‘^|| = 0{vN), the properties of 
the norm and Eq. (A25) imply 


= mn) , ( A 27 ) 

for n,m < cN. Additionally, conjugation by the corre¬ 
sponding Fourier transforms gives 

l(^™l(/yi^n)-(^m|p'|V'n)l = 

= \{f^i\{x^ym-{^|^m\£y^yn)\, (A28) 

and Lemma implies 

l(^ml(/)'l^n) - ■ (A29) 


Applying the triangle inequality to Eqs. (A27) and (A29) 
gives the desired result. □ 

Lemma 4. Given li,l 2 > 0, there exists a constant c > 0 
such that, for all n,m < cN, 

liAiiipA'Ax'^yyAt) - = i^an) . 

(A30) 

Proof. The property of the Hermite functions xtfnix) = 
y^{n+ l)/2tpri+1 (a::) + \An/2tf}n-1 ( x) immediately implies 

N/2-1 

x‘^\AA) = (27r/Af)i/‘^ X ^jAuixAlj) 

3 = -N/2 

= Vin+A/‘2\At+i) + vX^IV^n-i) > (A31) 

for all n > 0. Then, if x'-^ifnix) = Y.\A-i^G'An+i'{x), 
we obtain 

ix‘^y^\At)= X (A32) 

i'=-h 


(A24) and 


(AiiiP^yAxA^yAt) = X ^iAAi\iP‘^yyA^n+v) ■ 

i'=-h 

(A33) 

Because cu = a nd h and I 2 are constants, 

Lemma implies that Eq. (A331 can be approximated 
by 

I 2 

X ClAAm\p^yAn+l') = {Am\p^'-x’'^\An) (A34) 

V— — I 2 

within precision exponentially small in N, as long as 
n,m < cN. The constants c > 0 is as in Lemma [l] 
The constant /3 > 0 (used for the lower bound of i'i{N)) 
is as in Lemma [3l □ 
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A similar result is obtained if we swap the order of 
and p‘^, and x and p. This can be shown by acting with 
the corresponding Fourier transforms. 

Corollary 2. There exists a eonstant c > 0 sueh that, 
for all n < cN, 

||(7Jd-(n + l/2))|V-:J)f = z.i(fV). (A35) 

Proof. Alternatively, we can show that 

+ in + 1/2)^|^:J)| (A36) 

is exponentially small in N. Corollary implies that 
there exists c > 0 (as in Lemma[^ such that, if n < cN, 

|(V>::i(n+l/2f|V^(^)-(n+l/2f| =^i(7V) . (A37) 

Also, since + (p‘^)^)/2. Lemmas and 

imply {I = 2) 

\{'iPn\Hd\'fjn) - {i’n\H\tlJn)\ = | | i/d - (n + 1/2)| 

= MN), (A38) 

and then 

- in + 1/2)21 ^ ^ (^39) 


Since (iLd)^ = + (x‘^)2(p‘^)2 + 

(p‘^)2(a;‘^)2]/4, Lemma^(i = 2) implies 

mUHdrii’u) - in + 1/2)^I = i^diN) . (A40) 


It follows that Eq. (A361 can be approximated by (n + 
1/2)2 —2(n +1/2)2-|-(-n,-|-1/2)2 _ q precision that 

is exponentially small in N. 

□ 


Appendix B: High-order Trotter-Suzuki formula for 
the discrete QHO 

We first prove our results for the continuous variable 
QHO and then find approximations in the discrete case. 
Since [x,p] = i, the operators of the QHO form the Lie 
algebra sp{2) and satisfy the following commutation re¬ 
lations: 

[x‘^,f] = 2i{x,p} , [x^,{x,p}] = Aix^ , [f,{x,p}] = -Aif 

(HI) 


For s G R, Eqs. (HI) imply 


e p2g*s^ = p2 -p 2s{x,p} + 8s2i2 , 
^-^sp^^ 2 e^sp^ = _ 2s{x,p} + 8s2p2 , 

= {x,p} + Asx^ , 
g-*s33 {^^p}e*sp = {x,p} — Asp^ . 


(B2) 


We let U{t) = exp{—iHt) be the evolution operator 
of the QHO for time t G R. The second order Trotter- 
Suzuki (symmetric) approximation over a course of evo¬ 
lution time s is 

Ui{s) = ^ (B 3 ) 

While such an approximation is typically defined in a fi¬ 
nite dimensional Hilbert space, here we use it in Hilbert 
spaces of infinite dimension. We also construct higher or¬ 
der Trotter-Suzuki approximations using the recurrence 
relation 

C/p+i(s) = {Upisp)f Upis - 4sp) iUpisp)f , (B4) 

with Sp = s /(4 — 4 ^/( 2 p+i)) and p = 2,3,... [I5l ITHl 
IM| . The operator ep(s) = Up{s)U{—s) — 11, where 11 is 
the identity operation, can be used to quantify the error 
made in the approximation. 

Lemma 5. There exists a constant d > 1 such that, for 
all s satisfying |s| < 1/d and all n> 0, p > 1, 

||ep(s)|V'„)||=0((n + 2)|s|2p+i). (B5) 

Proof. With no loss of generality, we write 

ep(s) = [ ds' ds>ep{s')) . (B6) 

do 

Erom the definition of ep(s), we obtain 

a,.ep(s') = Up{s')fp{s')Ui-s') , (B7) 


where fp is an operator that depends on the approxi¬ 
mation order p. Furthermore, we can use Eqs. (B2) to 
show 

fpis) =ai(p)(s)x2-f6pp)(s)p2 + Q(p)(s){i,p} , (B8) 

where a/(p)(s), &/(p)(s), and c;(p)(s) are polynomials in s 
of lowest degree l{p), and l{p) is a positive integer that 
depends on p. For example, if p = 1, explicit calculation 
of Eq. (B7) results in 


+ e*""'/^(-fp2/2)e-**"'/^ - ix^/4 + iH , (B9) 


with P[ = {x^ -|-p2)/2. Equations (B2) imply /i(s) = 
—i(s‘^/4)a;2 — i(s^/ 2 )p 2 — i\s^/A){x'p\ and then l{p = 
1) = 2. The first goal is to obtain upper bounds on 
|ai(p)('®)l’ |ci(p)(s)|; we will obtain such 

bounds from the corresponding series expansions in s. 

In general, we will show by induction that l[p) = 2p. 
This result also follows from [T^. Since Up{s) = (1 -I- 


,{s))U[s), we can rewrite Eq. (B4) as 


Up+i{s) =((11 + ep(sp))H(sp))2(]l -p ep(s - 4sp)) 

H(s - 4sp)((]l + ep(sp))C/(sp))2 . (BIO) 
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This is a sum of 2® terms; the term without tp corre¬ 
sponds to {U{sp))^U{s — 4sp)(C/(sp))^ = C/(s), and the 
remaining terms sum up to C/p+i(s) — U(s). The sum of 
the terms containing a single Cp is 

Ep{s) =ep{sp)U{s) + U(sp)ep{sp)U{s - Sp)-|- 
-f U{^2sp)cp{s — 4sp)ff(s — 2sp)-t- 
+ U{s — 2sp)€p{sp)U {2sp) + 

+ U{s - Sp)ep{sp)U{sp) . (Bll) 


In the induction step we assume that l{p) = 2p, for 
some p > 1. A Taylor series expansion of ep(s) can 
be obtained by Taylor expanding Up{s') and [/(—s') in 
Eq. (B7). Because the lowest degree of /p(s') is assumed 
to be 2p, integration in s' implies that the lowest degree 
in the series expansion of ep(s) will be 2p-|-1. The lowest 
degree in the Taylor series of Ep{s) will be determined 
by the lowest degree of 4ep(sp) -|- ep(s — 4sp), which is the 
operator obtained from Ep(s) if we replace U hy 1 (i.e., 
the lowest degree term in the expansion of U). Then, 
since 4(sp)^P+^ + (s ~ dsp)^^''"^ = 0, the lowest degree in 
the series of Ep{s) is, at least, 2p + 2. Also, the lowest 
degree in the series of Up+i{s) — U{s) (and ep+i(s)) is de¬ 
termined by that of Ep{s) (i.e., the term of lowest order 
in Cp), and is also bounded from below by 2p + 2. It fol¬ 
lows that [/p_|_i(s) — [/(s) = -I-0(s^P+^), where V 

is some operator that depends on (powers of) p^ and 
{x,p}. The Trotter-Suzuki approximations determined 
by Eq. (B41 are symmetric and imply 


1 — Up+i{s)Up+i{~'S) 

= ([/(s) + + 0(s2P+3))x 

X ([/(-s) -b I>s2p+2 + o(s^P+^)) 

= 11 -b U(s)Vs^P+^ + Vs^P+^l7(-s) + 0(s^P+^) 

= 11 -b 2I/s2p+2 + 0(s2p-h3) ^ (3;^2) 


which can only be satisfied if E = 0. The last equality 
follows from the expansion of (7(s), whose lowest-degree 
term is 1. Then, the lowest degree in the Taylor series 
of [/p+i(s) - [/(s), or ep+i(s), is 2p -b 3 = 2{p -b 1) -b 1, 
implying that Zp+i = 2(p-b 1), and proving the induction 
step. _ 


From Eq. (B4), Up{s) is a product o f 0 (5^) exponen¬ 
tials of andp^ . The relations in Eqs. (|B2|) and Eq. (B7) 


imply that the highest degree in fp is smaller than 2 x . 
Since we already proved that l(p) = 2p, we obtain 


2x5*’ 

l—2p 


(B13) 


of p^ and {x,p}, needed to obtain fp, are linear com¬ 
binations of the same operators, and the largest prefac¬ 
tor in such combinations is a constant (8 in this case). 
The term of degree I in fp is obtained from 0{l) uni¬ 
tary transformations of the operators. Since |sp| < |s| 
and |s — 4sp| < |s|, there exists a constant d > 0 such 
that |tt;| < dK The value of d can be determined from 
Eqs. ( |B2[ ). If s is such that |s| < 1/d, the series for 
ai(p)(s) is convergent and |a/(p)(s)| = 0((|s|/d)^^’). That 
is, |a;(p)(s)| = 0(|sp^), and a similar result can be ob¬ 
tained for bpp){s) and cpp)(s). 

The properties of the Hermite functions imply 
ll^^l^n)ll = l|P^IV'n)|| = 0(n-b2), and also ||{x,p}|?/>„)|| = 
0{n + 2). From the triangle inequality and the previ¬ 
ous results, we obtain ||/p(s')|'*/'n)|| = 0{{n + 2)|s'p^'). 


In addition, sinc e U (—s')\ipn) 
im)('SOII = 1) Eq. (B7) implies 


e**'(”+i/2)|^/>„) and 


||ep(s)IV’n)|| < 



ds' ||/p(s')IV'n)ll=0((n + 2)|s|2^>+i), 

(B14) 


which is the desired result. 


□ 


Lemma ^ basically demonstrates that Up{s) is a prod¬ 
uct formula approximation of U (s) of order 2p -b 1 in s. 
The approximation is better for smaller values of n, i.e., 
for the low energy states. It is important to note that 
the dependence of the approximation error in n is only 
linear. 

Lemma 6. The number of exponentials of and p^ 
needed to prepare U{t)\'ipn), for |t| > 1, and within preci¬ 
sion e > 0, is 

M = e (|t| exp(7\/log((n-b2)|t|/e))j , (B15) 


where y > 0 is a constant. 

Proof. We first divide t into k intervals of size s = t/k, 
i.e., U[t) = {U{s))^. We will approximate each U{s) by 
Up{s), for some p > 1, and then choose p and k that mini¬ 
mize the number of exponentials in the product {Up{s))^ 
to obtain M. If |s| < 1/d, for some constant d > 0, 
subadditivity of errors and Lemma imply 

/fc((n + 2)|s|2P+i) = 0(e) . (B16) 


Each Up{s) can be implemented with less than 5^ expo¬ 
nentials of and p^. Then, the number of exponentials 


in {Up{s)p is Ad < kb^. Equation (B16l implies that 
there exists 


fc = 0 


(n-b2)i/2p|t|i-n/2p 


rl/2p 


(B17) 


and there is a similar expression for 6i(p)(s) and C;(p)(s). 
In order to show that the dominant term in ep(s) is that 
of degree 2p -b 1, as in the lemma, we still need to s how 
that the coefficients ui are bounded. In Eqs. ( |B2[ ) we 
showed that the corresponding unitary transformations 


to satisfy the desired error bound. This implies that the 
number of exponentials is, for any p = 1 , 2 ,..., 


M{p) = 0 


/5P(n-b2)i/2p|t|H-i/2p\ 


(B18) 
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It is simple to obtain the optimal value of p, defined by 
p = argminp/>i Ai(p'). The result is 

r\/log((n + 2)|t|/e)/(21og5)] >p, (B19) 

p > [v/log((n + 2 )|t|/e)/( 2 log5)J , (B20) 

implying that there is a constant 7 > 0 such that 

M = e (^|t| exp( 7 V'log((n + 2)|t|/e))) . (B21) 

The constant 7 can be obtained from the value of p and 
is approximately ^/2 log 5. The idea of computing an op¬ 
timal value of p was also considered in m- 

Note that |s| decreases with |t| so the assumption 
|s| < 1/d is valid, with no loss of generality. In par¬ 
ticular, we can always reduce s by a constant factor, at a 
constant increase in the cost, without changing the total 
order of operations. For example, we can assume that 
|s|5^ < c for any constant c = 0 ( 1 ), and still satisfy 

M = Q (||t| exp( 7 A/log((n-h 2 )|t|/e))). 

□ 


Three remarks are in order. First, so far we as¬ 
sumed h ^ 1 so we can disregard units. If the 
Hamiltonian is id = huj{x‘^ + p^)/2, and units are 
considered, then the number of exponentials is Ai = 


0 exp( 7 -\/log((n -I- 2)w|t|/e))^ . We also note that 
|s| 5 P = 0 ( 1 ) in Lemma Also, if |t| = 0(1) 
and e = 0 ( 1 ), the number of exponentials is Ai = 
0 (exp( 7 -\/log(n -|- 2))), for some constant 7 > 0. This 
implies lim„_>oo Ai/n“ = 0 and lim„_,.oo(logn)“/Ai = 0, 
for all a > 0. Finally, if s and p are as in Lemma[^ then 
||[(t/p(s))'' - [/(t)]|V'm)|| = 0 (e) for all m < n. 

In the discrete case, we define the N x N unitary ma¬ 
trices Up{s) by replacing x —>■ and p —» p‘^ i n th e 
dehnition of C/i(s) and Up{s) - see Eqs. ( |B3[ ) and (B4). 
The next goal is to use the previous results for contin¬ 
uous variables, and show that approximates 

(Op (s))*^|'!/'((), within precision 0 (e), for certain values 
of n. Because Up{s) is unitary, it will suffice to show 
(V'nKC^p (s))^l^n) is sufficiently close to (^/>„|([/p(s))*^|^/„) 
as these are close to 1 in absolute value. 


Lemma 7. Lets andp satisfy |s|5^ = 0(1). Then, there 
exists a constant c' > 0 such that, for all n < dN, 

\{i^n\UpisMi) - {llj„\Up{s)\lljn)\ = 0(5P(z/i(A))l"l) . 

(B22) 


Proof. By definition, Up{s) is a product of r < 5^ expo¬ 
nentials of {x'^Y {p'^Y- 

U^{s) = , ( b 23 ) 

The corresponding evolution times satisfy \ti\ < |s|. We 
first approximate each exponential by truncating the 
Taylor series at order I, where I will be chosen below. 


Since ||(x‘^)^|| = ||(p‘^)^|| = 0{N), the subadditivity of 
errors implies that the overall error in approximating 
all the exponentials in Eq. (B23) is 0{5P{N\s\y/l^-), as¬ 
suming that I > iV|s|. We then choose, for example, 
I = [2e|s|A] and Stirling’s approximation implies 

0{5P{N\s\y/ll) = 0(5P(1/2)2®I*I^) = 0(5P(j.i(A))l*l) . 

(B24) 


The property xYn{x) = \fnj2 '0n-i(a;) -I- 
\jynP l )/2 'ipn+iix) implies 

= VW^lYn-i) + Vin + l)/2|V’^+i) (B25) 

so that the approximation of ^ when acting on 

lYn)^ gives I'*/’«')• The coefficients c„/ can be 

obtained from the continuous-variable case; that is. 


n+2Z 


^{-itjx'^y /I'l IV'n) = ^ Cn'lYn') ■ (B26) 


Kl'=0 


n'—O 


li n + 21 < cN, for some constant c > 0 determined 
in Lemma and Lemma we can approximate \tpY) 
by lYn') within precision i^i{N) - see Eq. (A25). Since 
g-*L(p ) = pdg-itjix ) (^pyy, a truncated series ap¬ 

proximation of e“*b(p ) ^ when acting on I'i/')!), gives 
+ ^^i(A^)- The coefficients dn' can also 
be obtained from the corresponding continuous-variable 
case: 


C l \ n+2l 

'^{-itjfy'/l'\ I \Yn) = XI dn'\Yn') ■ (B27) 

V=0 / n'=0 

Because we approximate a product of 0(5^) expo¬ 
nentials, the above approximations are valid as long as 
n + 5^21 < cN. Then, the assumptions of working within 
the “low-energy” subspace apply in this analysis. Equiv¬ 
alently, we can assume that n < dN, for some constant 
d < c. The result is 

cN 

UyisMn) = X + Oib^IXiiNyP 

n'—O 

+0(5P(i.i(7V))l^l) . (B28) 

Because |s| = 0(1) and i 2 i{N) < 1, the dominant order 
in the approximation is 0(5 ^(j^i(A))I®I). The coefficients 
Xn' can be obtained from 

i i 

{'^{-itixyYh'-){'^{-it2fyyi2'-) ■ ■ ■ 

Li—0 I 2 —O 

I cN 

■ ■ ■ {'^i-drxy^/Ir'-Mn) = X ■ (B29) 

7=0 n'=0 

Since n + 5^(21) < cN, approximating Up{s) by trun¬ 
cating the Taylor series of each exponential of x^ and fd 
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at order I implies 

cN 

u^{s)Hn) = E Xr,,Wn) + o{bP{N\s\y/iy) 

n'—0 
cN 

= ^x„.K)+0(5"(^^iW)'^'). (B30) 

n'—O 

Then, 


Lemma 8. Let p be as in Eq. (B33), N' > 1, |t| > 1, 
and e > 0. Then, there exist eonstants c' > 0, 5 > 0, 
and dimension N = |’exp((5-\/log(A^'|t|/e)) + A^'/c'], sueh 
that 


\\[{U^{s)r-U^{t)Mu)\\=0{e), (B37) 

for all n < N', where k = t/s. The unitary {Up{s))’^ 
is a product of Ai = 0(|t|5^P) exponentials of and 

(pd)2. 


= {i^^\U^{s)\r^)+0{5^{MN)P) , (B31) 

where we used n < c'N < cN and |s| < 1, 5^ > 1. The 
result is 

(V'n|C/p(s)|V'n) = +0(5^K(iV))'*') • 

(B32) 

□ 


Proof p a nd s s atisfy 5^1 s| = c, for some constant c > 0 
- see Eq. (B33|. The subadditivity of errors, Lemma 
and Corollary imply 


l(^n|C^p(s)IV^n) - e' 


= 0{e\s/t\) + 0iN'\s\^P+^) 
= Oie\s/t\) , (B38) 


where we used Eq. (B16). This is valid if n < A^' < c'N, 


for some constant c' > 0, and our choice of N already 
satisfies such a condition. Then, using Corollary we 
obtain 


Eor the following results, we will assume that n < N' 
for some A^' > 1. We will set s and p as given by 
Lemma that is 

p = 0 (x/log(fV'|t|/e)) , (B33) 


and 5^|s| = c, for some constant c > 0 . 


Corollary 3. Let p be as in Eq. (B33), N' > 1, |t| > 1, 


and e > 0. Then, there e xist eonstant s E > 0, 5 > 0, and 
dimension N = |'exp((5-\/log(Af'|t|/e)) + N'/c] such that, 
for all n < N', 


\{'^n\Up{s)\ij^) - (V'„|C/p(s)|V'n)| = 0{e\s/t\) . (B34) 


IlC/p (s)IV'^) - = Oie\s/t\) + MN) . 

(B39) 

In particular, there exists a constant <5 > 0 such that, 
if A^ > exp((5-yiog(A^'|t|/e)), then vi{N) can be made 
negligible with respect to (e|s/t|). 

Using again the subadditivity property of errors, we 
obtain 

\m^{s)r\^u) - = 0(e) , (B40) 

for the corresponding choices of N and p. Also, Corol¬ 
lary [^implies, for n < N' < c'N, 

||(C/d(i) _ g-dn+i/ 2 )t)|^d^|| ^ o{\t\n,{N)N') . (B41) 


Proof. Note that n < N' < c'N, so that Eq. ( |B32 ) is 
valid. We then consider the term 0(5 p(i^i(A^))I^I)) in 
Lemma To make this term 0(e|s/t|), it sufhces to 
satisfy 


5Pg-/3Ar|d = 0(e|s/t|) = 0(e5-7|t|) , (B35) 

for some constant /? > 0 since vi{N) is exponentially 
small in N. This implies = 0(e/|t|) or, equiv¬ 

alently. 


/3Af|s| — 2plog5 = n(log(|t|/e)) . (B36) 

Since p^ = 0(log(A^'|t|/c))j ^ choice of N that satisfies 
/3A^|s| — 2plog5 = U(p^) also works. This can be sat¬ 
isfied if A^ = 0{p^/\s\), which is exponentially large in 
p. Then, t here exists a constant 5 > 0 such that, if 
N < exp((5-\/k)g(W]ti7^), the desired error bound is ob¬ 
tained if, in addition, N' < c'N. □ 


Then, there exists a constant i5 > 0 such that, if N > 
exp(5A/k)g(W]t|7e)), then 0{\t\L'i{N)N) = 0{e). The 
triangle inequality gives 

\\[{U^is)r-U^m^n)\\=0{e) (B42) 

for the corresponding choices of p and N, which deter¬ 
mines the first result of the Theorem. 

By definition, Up{s) contains less than 5^ exponentials 
of and The total number of exponentials 

in {Up{s))^ is bounded by \t/s\LP = 0{\t\5‘^P), for our 
choices of s and p. That is, there exists a constant 7 > 
0 such that the total number of exponentials is A 1 = 
0{\t\exp{-fy/\og{N'\t\/e))). 

□ 


Appendix C: Error bounds for the preparation of 
quantum states with Gaussian amplitudes 


We can use the previous Lemmas to prove one of our 
main results: 


In this section we prove some results regarding 
the preparation of quantum states with approximate 
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Gaussian-like amplitudes that serve as a basis for prepar¬ 
ing approximations of other eigenstates of the discrete 
QHO. The results in this section may be of independent 
interest. 

Lemma 9. Given N and 5 > 0, let 

N/2 

|y,d) = ( 1 / 0 ^) ^ exp(-jV(2(5)) \j) , (Cl) 

j=-Nl2 

with K = = Stt/N, t = 

\/(t^(2 — 4(t^)/2, and t' = l/(4t-|-dcr'^/t). Then, there 
exists a constant A > 0 such that 


^-ia(t)^ix t ^ip show that a similar result is ob¬ 

tained, up to some bounded approximation errors, when 
we work in the discrete Hilbert space of dimension N. 
We rewrite 

\^) = |(/j) , 

(C 8 ) 

where F is the operator that implements the Fourier 
transform, and we used the property * |^o) = 

g-ix t algo define the functions 


g{x) = (a;| F\(p) = ( 2 (T^/ 7 r)^/^ exp(— 


(C9) 


lll^d) _ g*a(t)e*G ) t e^P ) ^ o(exp(-A5)) . 


(C2) 


The global phase a{t) € R only depends on t. 


Proof. First we show that in CVs, Eq. (C2 1 holds exactly 


if the discrete states and operators are replaced by their 
CV versions. Then, we approximate the integrals ap¬ 
pearing in CVs by finite sums that appear in the actual 

Eq. ICl. 


where we used tpoix) = e~^^. Then, Eq. (C 8 ) 
becomes 

f dx /i(a;)e“ *g{x) = ^ — [ dx exp{—ax^) , 

J V’’■(I — *2t') J 

(Cll) 


In CVs, the wave function of a free particle evolves ac¬ 
cording to the Schrodinger equation, with a Hamiltonian 
—p^. If the initial state {t = 0) of the CV system is \ip) 
and the normalized wave function is 


and 


a = — it + 1/{2 — iAt') 

= [(7^ + 1/(2 + 8t'2)] - i[t - t'/[l + 4t'^)] . (C12) 


p{x,Q) = {x\ip) = 


(27ra2)i/4 

then, the wave function at time t > 0 is given by 


exp(—a:^/(4cr^)) , (C3) Next we approximate Eq. (Cll) by a finite sum, as- 


, . , 2cr^ 

p{x,t) = - 


2\ 1/4 


1 


\/—i2t + 2a‘^ 


exp{—x‘^/{—iAt + Aa‘^)) . 

(C4) 


We note that |(/ 3 (a;, 0 )p and \(p{x,t)\‘^ correspond to nor¬ 
mal distributions of variances cr^ and +t'^/a'^, respec¬ 
tively. We can rewrite (p{x,t) as 


suming the discretization where Xj = jA^/^ttT/V, j = 
—N/2,...,N/2 — 1, and bound the errors of the ap¬ 
proximation. Later, we relate this approximation with 

Eq. dC^ . 

We will show that one of the dominant sources of error 
in the approximation, for V 3> 1 , results from bounding 
j so that |j| < N/2. To estimate this approximation 
error, we consider 


Si — 




{2TT{a‘^ +t‘^/a'^)YN 

where the phase factor satisfies 


exp>{—x^/{Aa'^+At^/a‘^)), (C5) 


(2a2)i/4 ^ 


(C13) 


7 (a;,t) 


(C 6 ) 


In particular, 

OO 

I E 

j=N/2 


^ E e 

j=Nl2 


-U{a)nj ^ 


-5R(a)7riV/2 
1 _ g—SR(a)7r 


(C14) 


and a{t) solely depe nds on t, and can be computed ex¬ 
actly from Eq. (C4|. Selecting cr^ = Stt/N, for some 
given N > 0, and t so that 4 ct^ -|- = 2, we obtain 


where 3fi(a) is the real (and positive) part of a, and we 
used Xj > Xjqj^ = \/ttN/2. Then, since 3fi(a) > cr^ = 
Stt/N , 


if-1 


Ul e " e"'' " \g3) = (a:|' 0 o) = V'o(a;) 


(C7) 


with t' = l/(4t -I- Aa'^/t). Then, in CVs, we can ex¬ 
actly prepare the ground state of the CV QHO from 
the initial state |(p) by applying the unitary sequence 


I E 

j=Nl2 


G\<2- 


„-S^G2 


(CIS) 


(Stt^/N) ’ 

where we used that 1 — e~^ > x/2 ioi x < 1, which in 
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this case requires N > . We also note 






< 


(2a2)V4 


^/irt'N 

< 2{a'^/{TTN)f/'^ . (C16) 

Combining this with Eq. ( |C15 ), we obtain 

£1 < = 0(exp(-A(i)(5)) , (C17) 

where we replaced = 1 : 6 /N and > 0 is a const ant. 

We consider now the approximation of Eq. (Cll) by 
the infinite sum 


(2ct^)^/'^ ^27r ^ ^ 

v/^(l - i2t') V ^ Y 


(CIS) 


the last condition implies |(^'^) = b)j 

with (p{x) = + kT,0), i.e., the convolution of 

ip with the Dirac comb. Equivalently, we are defining 

N/ 2-1 

l<^")= E (C23) 

0=-N/2 

as used in Corollary Appendix]^ We note that 
p{x) - p{x,0) = E p{x + kT, 0) 

kpo 

= 0{N^/^exp{-TTN/{2a‘^))) (C24) 

and then there is a constant A*-^^ > 0 such that 

^3 = |||<^'') - E ‘^(^4.0) \j) II = 0{N^/^ exp(-A(3)iV2)) . 

3 

(C25) 


Using similar tools as in previous proofs based on the 
Dirac comb of period T = \/2nN, we obtain the error of 
the approximation of the integral as 


£2 = 


(2a2)i/4 


■\/2TTa{l — i2t') 


E' 

fc/O 


-uilKia) 


(C19) 


This error was obtained using the Fourier transform of 
exp(—az^). The frequencies are Wk = k'/^nN. With our 
choices for a, t, and t', we can bound 

S2 < (g^/(7r|a|))^/^ E 




^E' 

k^O 

<Y^e-N^k^ 

k^O 


(C20) 


where we first used |a| > 3fi(a) > and then 3fi(a) < 
3g^, which is valid under the assumption that N > 
or < I/tt. Then, there is a constant A^^^ > 0 such 
that £2 = 0 (exp(—A(^^iV^/i5)). 

The next goal is to show that the finite sum 


j=-Nl2 




3) ’ 


(C21) 


which we showed is a good approximation to Eq. (C 81 , 


also approximates the desired inner product be tween dis¬ 
crete states appearing in Eq. (C2). Equation (C21) can 
be realized as the inner product 

_ (C22) 

with (1‘^lF^'^lj) = h{xj) and (j|F^^|(^‘^) = g{xj). Using 
the properties of the discrete and CV Fourier transforms. 


We note that J2j \j) = (-^/(27r^(5))^/^v^|i^‘^), as 

defined in the statement of the Lemma. Similarly, the 
condition on h{xj) implies 

7V/2-1 Ar/2-1 

11'")= E ij)(Qie-“'‘'iV'o)= E 

j=-N/2 3 = -N/2 

(C26) 


with 




1/4 


(C27) 


being the convolution with the Dirac comb. If ^(a;) = 

||(x) - a^)\ = 0(exp(-D(7V))) . (C28) 

This implies that there is a constant A^"^^ > 0 such that 

£4 = IIIC'^) - E^(^/) !■?’) II = 0{Nexp{-\^*'>N)) . 

3 

(C29) 

We then define |^'^) = ^i^j) \j) using the defini¬ 

tion of IV'o), we obtain |^'^) = (jV/(27r))^/'^e~d^‘^) ^* |'0 n)- 
In summ ary, so far we demonstrated that Eq. ( |C 8 | ), or 
Eq. (Cll), can be approximated by 


/I \ 

[-s) (C30) 

within approximation error of order £i-|-£ 2 -l-£ 3 -l-£ 4 , which 
can be obtained by using the subadditivity property of 
errors. 

There are two additional approximations that still need 
to be bounded. One is £5 = IH'i/'o) ~ (-f"^)^|i/’o)ll- Using 
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the results in Eq. (A25), we can show that |^/>q) is al¬ 
most an eigenstate of of eigenvalue. In particular, 
Eq. (A25) implies £5 = 0(exp(—/3A^)), for some constant 
/3 > 0. Last, we seek to show that k approximates 
We write VttS = \JNj^-K f dx which can 

be approximated by the infinite sum e~^ within 
approximation error 


£6 = v^v^E' 












(C31) 


where Wfc = fc-\/27rA^, and Eq was obtained by comput¬ 
ing the Fourier transform of a Gaussian with variance 
ttS/N. Then, there is a constant > 0 such that 
Eq = 0(5exp(—A^®^i5)). Setting a cutoff in the infinite 
sum so that \j\ < N/2, this introduces an additional 
approximation error £7 = 0 (exp(—A^^^A^^/5)), for some 
constant A^^^ > 0 . 


The overall result is that we can approximate Eq. (C 81 , 


or Eq. (jC^, by 

(4^1 , (C32) 

within precision But since = 

0 , and F^)'^ = I, the above equation is exactly 

(' 0 o|e®*'® ^ ^ *|E)) desired quantity for the 

Lemma. In particular, this implies 




- (V'i 




IE 




"t'gdpOA 


lE 


<E^* 


i=l 


(C33) 

(C34) 

□ 


Appendix D: The discrete Jaynes-Cummings 
Hamiltonian 

The discrete Jaynes-Cummings Hamiltonian was 
introduced in Eq. (38). The Hilbert space dimension 
is 2N, where N is the dimension for x^ and p®. In 
this Appx., we first prove that the normalized states 
l7l±) = E[IE)|0) ± IE-hi)|1)]: where n < cN for 
some constant c > 0 , are almost eigenstates of of 
eigenvalues ±y/n -I- 1 . 

Lemma 10. There exists a constant c > 0 such that 

\\Hfc - (±E^) 17^)11 = ^i(^) > (Dl) 

for all n < N' < cN, where vi{N) = exp(—H(A^)). 

Proof. First, we note |0) = (a;'^ — |I). 

As in Eq. (A3I), satisfies, for all n< N, 


= EE + i)/2|E-n) + EE2 |E-i) ■ (D2) 


Lemma implies 

llElE) - EE + i)/2|E-n) + *EE2 |E-i)II = 

(D3) 

for all n < N' < cN, where 1 > c > 0, and i'i{N) = 
exp(—f2(iV)). Then, 


IliLElE) | 0 ) - E^lE+i) |1) II = MN) (D4) 
Ili^ElE+i) |i) - E^IE) |o) II = , (d5) 

where the last equation follows by doing a similar analy¬ 
sis. Since ||H^cll =0{N^^‘^), we can replace by the 
actual eigenstates |(()®) using Eq. (19), and the order of 


the approximation error is also exponentially small in N. 
It follows that 


IIhE - (±E^)|7(!,±)II = ^i{N) . (D 6 ) 

□ 

Next, we seek efficient decompositions of the evolu¬ 
tion operator induced by For simplicity, we use an 

asymmetric first order Trotter-Suzuki approximation as 
this already provides the desired scaling with TV, how¬ 
ever, one may use higher-order approximations to im¬ 
prove upon the scaling on the precision parameter. 


Lemma 11. Let 

W{s) = e 


— p-i(2:‘*®o’i:)s/-v/2gi(p‘*0o-„)s/-v/2 


(D7) 


Then, there exists a constant c > such that, given n < 
N' < cN and s = Oiejsjn F 1), 

||(H/(s)-e-^-^)|E)|0)||=O(e^). (D 8 ) 

Proof. We first define £i(s) = IT(s)e®^Jc® — il, and, since 
£i(s) = ds'ds'Ei{s'), we obtain 


£i(s)= I ds'IT(s')E(^) 


/o 


k>l 


X [(p® (g) ay), [(p® ® ay), [..., (g) ay)].. . 

(D9) 

Because ||a:‘^|| = ||p®|| = 0{'/N), we can cutoff the sum 
in fc at AT = Of\/N) with an approximation error that is 
negligible (i.e., exponentially small in ^fN). Since K < 
cN, for some constant c > 0, we then use Eqs. (D2) 
and (D3) to show 


= 0(E(n -I- 1 )... (n -I- k)) , (DIO) 

where = k < K. Since |7®±) approximate 

eigenstates of when n < N' < cN, and these are 
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combinations of |(/)^) |0) and \(j)n+i) |1), we obtain 

Iki(5)l7;!.±)ll = 


If s < ejy^2(n + 1), Stirling’s approximation implies 


= o( [ ds y (( s '^/ 2 ) Vfc!) + 2 ). .. (n + A: + 1) I /k \ 

J lki(s)|7;j,±)ll = O , (D13) 

/ K \ \k=l / 


= O ( J2{s^+\V2)^/{k + l)!)v/(n + 2)...(n + fc + l)j 

(Dll) 


^fc=i 


where we assumed that other contributions to the error 
that are exponentially small in or VlV are negligible in 
the asymptotic limit. We use the bound on the binomial 

coefficient to obtain 


and then ||£i(s)|7^_±)|| = O (e^). 


□ 


K 


ki(s)|7)!.±)|| = O 1 + k + 

(D12) 





